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Series Foreword 


Dear reader, 

the series Simula SpringerBriefs on Computing was established in 2016, with 
the aim of publishing compact introductions and state-of-the-art overviews of 
select fields in computing. Research is increasingly interdisciplinary, and students 
and experienced researchers both often face the need to learn the foundations, 
tools, and methods of a new field. This process can be demanding, and typically 
involves extensive reading of multidisciplinary publications with different notations, 
terminologies and styles of presentation. The briefs in this series are meant to ease the 
process by explaining important concepts and theories in a specific interdisciplinary 
field without assuming extensive disciplinary knowledge and by outlining open 
research challenges and posing critical questions in the field. 

Simula has a major research program in computational physiology that includes 
a long and close collaboration with the University of California (UC) San Diego. To 
reflect this research focus, we established in 2020 a new subseries entitled Simula 
Springer Briefs on Computing - Reports on Computational Physiology. The subseries 
includes both introductory and advanced texts on select fields of computational 
physiology, designed to advance interdisciplinary scientific literacy and promote 
effective communication and collaboration in the field. This subseries is also the 
outlet for collections of reports from the annual Summer School in Computational 
Physiology, organized by Simula, University of Oslo, and UC San Diego. The school 
starts in June each year with students spending two weeks in Oslo learning the 
principles underlying mathematical models commonly used in studying the heart 
and the brain. During their stay in Oslo, students are assigned a research project to 
work on over the summer. In August, they travel to San Diego for another week of 
training and project work, and a final presentation of their findings. Every year, we 
have been impressed by the students’ creativity and we often see results that could 
lead to a scientific publication. Starting with the 2021 edition of the summer school, 
we have taken the course one step further by having each team conclude their project 
with a scientific report that can pass rigorous peer review as a publication in this 
subseries. 
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All items in the main series and the subseries are published within the 
SpringerOpen framework, as this will allow authors to use the series to publish 
an initial version of their manuscript that could subsequently evolve into a full-scale 
book on a broader theme. Since the briefs are freely available online, the authors do 
not receive any direct income from the sales; however, remuneration is provided for 
every completed manuscript. Briefs are written on the basis of an invitation from a 
member of the editorial board. 

Suggestions for possible topics are most welcome, and interested authors are 
encouraged to contact a member of the editorial board. 
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Preface 


Why would you want to read these notes? 


When something is very large, very small or very complex, it's often helpful to 
represent it using a model in order to understand what is going on. Nature offers 
many examples of situations where direct interrogations are difficult, and a model 
can provide insight into the phenomena of interest. Physiology, for example, is a 
field that is complex and sometimes poorly understood, and where mathematical 
models are frequently used to increase understanding. These models are often 
formulated in terms of differential equations. However, many students who want 
to learn state-of-the-art physiology may not be familiar with differential equations. 
They are probably able to use software tools that solve these equations in numerical 
simulations, but they may not fully understand the underlying principles. Our goal 
with these notes is to provide a simple introduction to differential equations and 
give examples of how to solve them. Differential equations is a vast area of active 
research, and we can only provide a glimpse. But, we hope that these notes will 
provide a basic understanding of the principles, such that computational codes will 
appear less like a black box, and more like something that is comprehensible. 


The Simula Summer School in Computational Physiology 


Every year Simula Research Laboratory organizes a summer school in Computational 
Physiology together with the University of California, San Diego (UCSD) and the 
University of Oslo (UiO). The students come to Oslo in June and learn about models 
and software used in computational physiology. Next, the students are divided into 
groups and work on different projects through the summer. АП groups are assigned 
one or more mentors who help the students with their investigations. In August, all 
students and mentors meet again at UCSD and continue their work on the project. In 


x Preface 


addition, they attend guest lectures at UCSD and a two-day workshop on scientific 
writing organized by experienced editors of Nature. 

The models and software used in the course are state-of-the art and the students 
complete the course by writing a scientific report that is published in a Simula 
SpringerBriefs on Computing in the sub-series on Computational Physiology. In 
order to reach the advanced level of using state-of-the-art methods and software, there 
is not enough time to cover the details of all the subjects covered in the course. The 
students are most often enrolled in MSc., PhD or Post Doc programs in universities 
around the world. Their scientific backgrounds vary greatly. There are students 
from theoretical disciplines like computer science, scientific computing, statistics or 
mathematics, and from more classical science disciplines like physics, biology, or 
medicine. And perhaps the largest group of students comes with a background in 
bioengineering. We notice that the students easily follow 'their' part of the course 
and struggle with the parts they are unfamiliar with. 

Scientific work in computational physiology is inherently interdisciplinary so 
meeting this reality in the summer school is proper training. Typical research projects 
comprise elements of many disciplines and thus it is very common to not fully 
understand all elements of a project, and almost no project is completed by only 
one person. Similarly, it is exceptionally rare to see a single-authored paper in 
computational physiology, or in computational science in general for that matter. So, 
learning to communicate across disciplines is very useful, but also very challenging. 

The most common mathematical machinery used to model physiology (and 
physics in general) is differential equations. Most of the models introduced in 
the summer school is founded, in one way or another, on differential equations. 
Elements of this subject are taught at every university in the world and students 
from mathematics, physics or scientific computing most likely have a course in this 
subject. But students from biology, medicine or even computer science rarely know 
much about differential equations. Now, the summer school is very streamlined and 
the software is usually prepared and can be used without a complete understanding 
of the underlying models, but it is clearly advantageous to know, at least intuitively, 
the foundation of the models and solution methods that are used in the course. 


How are the notes organized? 


We have organized these notes in two parts. In the first part (Chapters 1—7), we 
introduce the concept of differential equations and numerical methods. To keep the 
exposition as simple as possible, we use simple, unitless equations as examples 
in these chapters. In the second part (Chapters 8-12), we apply the techniques 
introduced in the first part of the book to a selection of models of electrophysiology. 

More specifically, we start these notes by considering a very simple differential 
equation. For this equation, we introduce a numerical method and we study the error 
of the method. The simplicity of the equation allows us to study these important 
concepts in a very explicit manner. Next, we move on to systems of ordinary 
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differential equations and we show how the FitzHugh-Nagumo model can be solved 
numerically. Partial differential equations are then introduced and we show how the 
diffusion equation can be solved numerically. By combining the diffusion equation 
and the FitzHugh-Nagumo model, we introduce traveling wave solutions and show 
how a reaction-diffusion system can be handled numerically. 

At the beginning of the second part of the notes, in Chapter 8, we introduce 
the Hodgkin-Huxley equations and from this chapter we start using units for all 
quantities involved. The Hodgkin-Huxley equations are the most famous system of 
equations in physiology. The equations model the action potential of an axon and 
have been proved to represent that process with great accuracy. We also consider 
a similar model for a cardiac action potential. The membrane models of cardiac 
electrophysiology have evolved into a very complex matter, but the structure remains 
very similar to the Hodgkin-Huxley model. 

After being familiarized with the Hodgkin-Huxley membrane model, we combine 
it with the cable equation to model the propagation of an action potential along an 
axon. Next, we introduce the bidomain and the monodomain models in two spatial 
dimensions. The monodomain equation is very similar to the cable equation, but the 
numerical solution will be a bit more complicated because we consider a 2D problem. 
The bidomain model is more complex than the monodomain model, but we will see 
that by using operator splitting, solution methods become quite straightforward. 
Never heard of operator splitting? You will learn about it from these notes. It is 
a very powerful technique used to break complicated problems into problems we 
already know how to deal with. We will explain it in Chapter 7 and show some 
applications later in the notes. 

One major issue with the monodomain and bidomain models is that they represent 
averaged quantities in the sense that the cardiomyocytes are not present in the model. 
In fact, both the extracellular (E) space, the cell membrane (M) and the intracellular 
(I) space are assumed to exist everywhere. That is a bold assumption, so we will also 
introduce the EMI model where all these elements (E, M, and I) are explicitly part 
of the model. 

The monodomain and bidomain models provide reasonable approximations 
of cardiac electrophysiology at the mm-scale, and the EMI model addresses 
electrophysiology at the wm-scale. The next level is the nm scale!. The relevant 
equations at the nm-scale are the Poisson-Nernst-Planck (PNP) equations. From a 
computational perspective, the PNP equations are extremely challenging and can, at 
present, only be used to study very small regions; one complete cardiomyocyte is 
far too large to be modeled by the PNP equations. We will show how to solve these 
equations in the final chapter using operator splitting and finite differences. 

In order to keep these notes relatively brief, many details about the methods 
and models introduced must inevitably be left out. At the end of most chapters, we 
therefore include a section called 'Comments and Further Reading’, providing some 
comments on these details and suggestions for further reading. 


! Note here that meter is denoted by m, and mm means millimeter or 107?m, иш is micrometer 
and means 107m, and, finally, nm means nanometer or 107?m. If you ever need help with units, 
we recommend www .wolframalpha.com. 
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Why do we use the finite difference method? 


There are many ways to solve differential equations. In computational physiology, the 
finite element method may very well be the most popular alternative. The reason for 
this is that the geometry ofthe computational domain (e.g., the heart) is quite complex 
and thus very difficult to represent using a finite difference method. Finite difference 
methods are easiest to work with when the computational mesh is very regular, 
but the finite element method is constructed to allow for highly irregular meshes. 
Both finite element methods and finite volume methods are successfully applied to 
simulate complex phenomenas on complex geometries. The finite difference method 
is successfully applied to simulate complex dynamics but the code quickly becomes 
clunky in the presence of geometries that are any more complex than a cuboid. 
Nevertheless, we have chosen to focus entirely on the finite difference method in these 
notes and the reason for this is simplicity. The method is more or less completely 
defined by simply replacing derivatives by differences. Therefore, we stick to simple 
geometries and use finite differences. In the summer school, much more advanced 
simulations, using finite elements, will be performed, but we still think it is useful 
to know how this in principle can be done using the simplest possible method. 


It's open access 


These notes are printed in Simula SpringerBriefs on Computing in the sub-series on 
Computational Physiology. That means that the notes can be downloaded for free. 
The software (Matlab) used to generate all figures and tables in these notes are freely 
available online. The codes are written primarily for clarity and less emphasis is put 
on efficiency. If you find a bug in the codes or an error in these notes, please send a 
mail to aslak 9 simula.no. The software associated with these notes can be found at: 
https://github.com/karolihj/differential-equations-book- 2023 
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Part I 
Tools for Differential Equations and 
Numerical Methods 


A 
Check for 


Chapter 1 weirs 
Getting Started 


In this introductory chapter we will introduce two essential concepts: differential 
equations and numerical methods. If you want to spend as little energy as possible 
(don’t be ashamed of that - energy preservation is both fashionable and a fundamental 
property of many biological mechanisms - it’s fine) you can get a good overview 
from this chapter alone. 


1.1 What Is a Differential Equation? 


You have no doubt seen an algebraic equation. A typical algebraic equation may 
look like 
x?-4x 4320 (1.1) 


with solutions x = 1 and x = 3. So, solutions of algebraic equations are numbers. 
The solutions of differential equations, on the other hand, are functions. One very 
simple differential equation is given by 


y'(t) = y(t). (1.2) 


Here, we typically assume that t represents time, and that the equation (1.2) describes 
how the function y changes with time. Suppose we also know that the solution is 1 
at f — O, that is 

y(0) = 1. (1.3) 


Such a condition is generally referred to as an initial condition and is needed in order 
to find a unique solution of the problem. In this case, the solution is given by the 
function 

y(t) = е’. (1.4) 


It is straightforward to verify this. We simply note that у(0) = e°? = 1 and, by 
differentiation, that y’(t) = e' = y(t). So all is good; we have found the solution of 
our first differential equation. 
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For a long time, differential equations were solved in this way. Formulas were 
derived using pencil and paper. If the equations were very complex, they were 
simplified in order to be solved by a formula and people spent whole academic 
careers deriving such approximate solutions of differential equations. The reason 
for this was that the solutions of the equations could bring critical insight into a 
phenomenon modeled by the equation. Nowadays, differential equations are solved 
by computers. In some cases analytical formulas can be derived and then it is usually 
done by computing systems like Mathematica or Maple, or other tools for symbolic 
computations. But the more common approach is to solve the equations numerically. 
In order to do that, we need to transform the equations into a form that is suitable 
for computers. Almost all differential equations that are solved in computational 
physiology are solved by computers. It is the main purpose of these notes to teach 
you the basics of how that is done, and we might as well start with the very simple 
case that we've already introduced. 


1.2 What Is a Numerical Method? 


A numerical method is a way to solve a mathematical problem on a computer. We 
will see that one way to prepare differential equations for solution on computers is 
to replace derivatives by differences. That may not come as a surprise since you may 
remember from calculus that the derivative was introduced as the limit of a difference. 
Concretely, for a smooth! function f = f(t), the derivative is, per definition, 


/@ hn - fo) 


A: (1.5) 


'(t) = li 
ше да, 
Instead of going all the way to zero, we can settle for a small value of At and use the 


approximation 
pays A507 10. 
At 
Here, t is time, and we are interested in solving the equation from t = 0 (where the 
solution, y, is known to be 1) until some t = T, for example T = 1. For some fixed 
value of At, it is useful to define a number of discrete points in time, 


(1.6) 


th = n X At, (1.7) 


for n = 0,..., N, where N = T/At. We note here that the step in time from f, to tn 
is At, and that the final time is given by ty = N x At = T. We want to compute a 
numerical approximation to the solution of the problem 


! A lot could be said about smooth functions. Generally, regularity of functions is important in 
solving differential equations. But for these notes it is sufficient to simply think of smooth functions 
as — smooth. 
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y(t) = y(t), (1.8) 
y(0) = 1, (1.9) 


at the time steps {t, |4 о by replacing the derivative of y by the formula given in 
(1.6). This replacement can be written as 


Уп+1 = Yn 


= Уп, 1.1 
dh ly (1.10) 


уо = 1, (1.11) 


where, in general, у, denotes an approximation of the solution of (1.8) at time tn; 
i.e., ул = y(t,). By rearranging (1.10), we find that 


Уп+1 = (1 + At)yn, (1.12) 


and we refer to this as the computational form of the numerical scheme. Since yo = 1, 
we find y; = 1 + At, and y2 = (1 + At)y; = (1 + Ap)’, and so forth. In general, we 
have 


уп = (1 At). (1.13) 


So in this unusually simple case, we don't need a program to compute the numerical 
solution. That is very unusual and we use this example just to introduce the concepts. 
More commonly, we first need to compute the solution y; from the initial condition, 
yo, using a formula like (1.12), and then insert this solution y, into (1.12) to find 
the solution y2, and so on until we find yy using the previously computed ум. 
This type of repetitive computation is usually most conveniently performed using a 
computer program. 

In Fig. 1.1 we have plotted the analytical (i.e. exact) solution, y, together with the 
numerical solution, y,, for some different choices of the time step, At. We observe 
that when Aż is small, the numerical and analytical solutions are indistinguishable, 
but that for larger values of Ат, for example At = 0.2, there is a visible difference 
between the two solutions. The analytical solution is exact, which means that the 
difference between the two solutions are due to an error of the numerical method. In 
the next section, we will consider this error in some more detail. 


At = 0.01 3 At = 0.001 


analytical 
2.5 2.5| MS numerical 


0 0.5 1 0 0.5 1 


Fig. 1.1 Analytical (solid line) and numerical (dotted line) solutions of the differential equation 
(1.8)-(1.9) for t between 0 and 1 for some different values of Ar. 
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1.3 What Is the Error of a Numerical Method? 


The error of a numerical method is the difference between the numerical and the 
exact solution of the problem. Since we are in the fortunate position of knowing 
both the analytical (exact) solution, y(t), and the numerical solution, y„, we can also 
compute the error? of the numerical solution defined by 


En = |y(tn) – уһ. (1.14) 


As observed in Fig. 1.1, the error depends on the size of the time step, At, which 
is natural since we are supposed to pass to the limit of Ar = 0 in the definition of 
the derivative (see (1.6)). When we approximate 0 by something small, we must be 
prepared for it to come with a price — and the price is the error we will encounter. To 
illustrate this, we put 7 = 1 and compare the analytical solution y(T) = e with the 
numerical solution at 7 = 1 for several values of At. Since the final time is T = 1, 
we must have N x At = 1 and therefore we compare 


1 N 
"T (i + А) (115) 


with the exact solution y(1) = e. In Table 1.1, we show the error Ey = |y(1) – умі 
for N =5, 10, 100, and 1000. We also show En /At and we note that this value seems 
to be more or less constant. It is thus evident that when Af goes to zero (N goes to 
infinity), the numerical solution converges towards the analytical solution?. 


Table 1.1 Error of the numerical solution of the differential equation (1.2)-(1.4) at t = T = 1 for 


different values of At — A. The error is defined as Ey = |y(1) – yw |, where y(1) = e is the 


analytical solution, and yy is the numerical solution, defined by (1.15). 


N At EN En /At 
5 0.2 0.23 1.15 
10 0.1 0.125 1.25 
100 0.01 0.0135 1.35 
1000 0.001 0.00136 1.36 


? 'The error defined here is referred to as the absolute error. An alternative is to consider the relative 
error given by 
ly (tn) — Yn | 
bu) ^ 


3 [f this feels familiar, it may be because you learned in calculus that 


lim(1 + е) = е. 
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We have learned that by replacing the derivative y’(t) by a finite difference 
approximation (yn+1 = Yn)/At we can find an approximate solution. Usually, the 
finite difference scheme must be implemented on a computer but in this very simple 
case, we can find a formula for both the numerical and the analytical solutions. 
The error introduced by this method is Ey = 1.36 x At. This indicates a relation 
that is generally true: The error is smaller for smaller values of Ат, but the work 
associated with running through the time steps to compute the numerical solutions 
is increasing for smaller values of At. We will see plenty of examples of the fact 
that finer resolutions (i.e., smaller At) means higher accuracy and more work. So 
life is fair; a low quality solution is cheap and a high quality solution is expensive. 
Since Ey is proportional to At, we have linear (or first order) convergence. If the 
error had been proportional to Аг? we would have had quadratic (or second order) 
convergence, and so on. 


1.4 Implicit vs. Explicit Numerical Schemes 


At this point, we will briefly mention the difference between an implicit and an 
explicit numerical scheme. When a numerical scheme can be written in the generic 
form 

ул = Е(у"), 


like in (1.12), we refer to the scheme as an explicit scheme because y"*! can be 
explicitly computed as a function of y". Conversely, if an equation has to be solved 
in order to compute y"*! based on y", the scheme is referred to as implicit. Implicit 
schemes will be introduced in Chapter 4 of these notes. 


1.5 But What Is a Differential Equation? 


We have given you one example of a differential equation, and there are many more 
in the subsequent chapters. The easiest way to grasp what differential equations are 
is probably by seeing many examples. However, in general, differential equations 
are mathematical relations where the unknown is a function, and derivatives of the 
function are used in the formulation of the equation. The solution is a function. 
Classical mathematical questions related to differential equations are: Existence: Is 
there a solution of the equation? Uniqueness: Is there only one solution? Stability: 
Can we slightly perturb the parameters of the equation and obtain almost the same 
solution? And Solution: How can we find the solution? In these notes, our emphasis 
will be on the last question and we will concentrate on numerical solutions of the 
equations. 
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1.6 Analytical Solutions 


Our focus in these notes will be on numerical methods for solving differential 
equations. However, we will also encounter some cases where analytical solutions 
(i.e., solutions given by a formula) can be found. While finding analytical solutions 
to differential equations is a vast field, these techniques are not typically applicable 
to the types of problems we will be studying, which is why we will rely on numerical 
methods. In this section, we will indicate a couple of possible approaches to obtain 
analytical solutions to illustrate that it is not entirely mysterious, but we will not 
delve into the methods in these notes, and you can safely skip this subsection and 
jump to Chapter 2 if you want. 

First, we will demonstrate an approach for obtaining an analytical solution for the 
simple example studied earlier. We start by repeating that the differential equation is 
given by 

у'@) = у@) (1.16) 


with the initial condition y(0) = 1. To find an analytical solution for this problem, 
we first write the equation in the form 


y'(t) 
= (1.17) 
y(t) 
Now, we can integrate both sides from 0 to t and get 
ty t 
li Aa | ldr. (1.18) 
o X(T) 0 
Therefore, 
[П(у(т))} = t. (1.19) 
ог 
In(y(t)) – In(y(0)) = t. (1.20) 
Since y(0) = 1, and thus In(y(0)) = 0, (1.20) reads 
In(y(t)) = t, (1.21) 


and raising e to the power of both sides of the equation, we get the analytical solution 
y(t) = et. (1.22) 


This way of finding the solution of a differential equation can be generalized to 
equations of the form 


F'(y)y'(t) = FOA) (1.23) 


with the initial condition y(0) = yo, where yo is a given number. Here, F = F(y) is 
assumed to be some invertible function. By following the steps above we find that 
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[In(F(y(7)))]o = t. (1.24) 
ү F(y(t)) 
y(t))\ _ 
| ПО) | ={ (1.25) 
SO 
FOH) = F(yo)e’, (1.26) 
and, finally, 
y(t) = F7'(F(yo)e"). (1.27) 


1.6.1 Integrating Factors 


Another way of obtaining analytical solutions of differential equations is by applying 
integrating factors. We can illustrate this technique by considering the equation 


y'(t) + p(t)y(t) = q(t), (1.28) 


with the initial condition y(0) = yo, where yo is given. Here, we assume that р = p(t) 
and q = q(t) are known functions. In addition, we assume that we have a function 
P = P(t) with the special property that 


P'(t) = p(t). (1.29) 
By multiplying (1.28) by the integrating factor 


eb) 


we get 
ePO y(t) + ePOp(r)y(r) = eP q(t). (1.30) 


Here we observe that (because of (1.29)) the left-hand side can be written in a more 
compact manner as, 


eP OYA)  eP?p(r)y(r) = (POYA. (131) 
and therefore (1.30) can be rewritten to read, 
(eP OOY = eF?q(). (1.32) 


Integrating both sides from 0 to г, we obtain, 


t 
eP O y(t) – eP y = ] еР©)а(т)ат. (1.33) 
0 


Finally, the solution is given by 
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t 

y(t) = e PO Gas +f "дат . (1.34) 
0 


In the case of the simple example equation considered in this chapter, y’(t) = y(t), 
we note that the equation can be written in the form (1.28) for p(t) = —1, q(t) = 0, 
and yo = 1. Furthermore, (1.29) is fulfilled if P(t) = —t. Inserting these functions 
into (1.34), we obtain the analytical solution 


t 
xo et [oie fen] =e (135) 
0 


1.7 Comments and Further Reading 


Here is a list of suggested further reading on a few of the overarching topics of these 
notes. 


1. Throughout these lecture notes, we will refer to results from calculus. A good 
book on the topic is [12]. 

2. Introductions to the basics of scientific computing are presented in, e.g., [3, 4, 
5,6, 10, 11, 145 17]. 

3. Introductions to differential equations can be found in [7, 8, 13, 16]. There are 
many texts on analytical solutions of differential equations; one comprehensive 
collection of solutions is given in [15]. Analytical techniques like the methods 
described above are well covered in [2]. 

4. A comprehensive introduction to mathematical biology can be found in [9]. 

5. For those interested in reading about the place of mathematics in biology, and 
the role of biology in mathematics, we strongly recommend [1]. 
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Chapter 2 e 
A System of Ordinary Differential Equations 


Above, we introduced a very simple differential equation given by 
y(t) = y(t). (2.1) 


This is referred to as an ordinary differential equation (ODE) because it involves 
the derivative with respect to only one variable. If the derivative with respect to 
more than one variable is involved, the equation is referred to as a partial differential 
equation (PDE). We will come back to PDEs in Chapter 3, and spend some time 
discussing how to solve them, but in this chapter, we will keep focusing on ODEs. 

The ODE (2.1) is a scalar equation because there is only a single unknown 
function to be found. Now, we will start considering systems of ODEs. A typical 
system of ODEs can be written in the form 


y'(t) = FOE). (2.2) 


Here, y is a vector and F is vector valued function. 


2.1 The FitzHugh-Nagumo Model 


We will introduce numerical methods for systems of ODEs by considering the 
celebrated! FitzHugh-Nagumo model published by FitzHugh [1] in 1961 and, 
independently, by Nagumo et. al. [2] in 1962. The model is a system of ordinary 
differential equations with two unknowns, and is commonly used as a simple model 
for the action potentials of excitable pacemaker cells. 

We consider the following version of the FitzHugh-Nagumo model, 


v’ = cqv(v — a)(1 — v) – cou, (2.3) 
w’ = b(v — du). (2.4) 


! These two papers together are cited more than 11,000 times. 
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Fig. 2.1 Numerical solutions v,, (left) and w, (right) of the FitzHugh-Nagumo model specified by 
(2.3)- (2.5). The numerical scheme used to compute the solutions is specified in (2.8)-(2.9), and 
we have used At = | and the initial conditions uy = 0.26 and wọ = 0. 


Here, the constants are given by 
a = —0.12, cı = 0.175, c2 = 0.03, b = 0.011, d = 0.55, (2.5) 


and the unknown functions are v and w. In order to solve the system of equations 
numerically, we use the steps introduced for the scalar equation in Chapter 1 and 
start by replacing derivatives by differences. The discrete system then reads 


Un] — Un 


At 
Unc] — Ш 


Xt = b(vn — dwn). (2.7) 


= C1Un(Un — a)(1 — vn) — C2Wn, (2.6) 


Again, we reorganize this system to write it in computational form, 


Un+1 = Un + Аі [с Un(Un а)(1 С Un) C2Wy], (2.8) 
Wn + At[b(v, — dwn). (2.9) 


Un 


This time, however, we note that we will need a piece of software to compute the 
solutions. But it is straightforward to implement this since the numerical solution at 
time tn+1 is an explicit function of the numerical solution at time ty. 


2.1.1 Numerical Computations 


We assume that the solution is known initially, so we define (for instance) vg = 0.26 
and wo = О. Here, it is useful to note that if we put both vp and wo equal to zero, 
the solution will remain zero (for both v and ш) for all time. But if we perturb v a 
little, we get very different solutions. In Fig. 2.1, we show the numerical solution 
from t = 0 tot = T = 5000. In the computation, we have used N = 5000 which 
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Fig. 2.2 The numerical 02 
solutions v, and w, of the 
FitzHugh-Nagumo model 0 
from Fig. 2.1 displayed in a 02 
parametric plot with vn on the 0.5 0 0.5 1 
x-axis and wn on the y-axis. v 


Table 2.1 Error of the numerical solution of the FitzHugh-Nagumo model specified by (2.3)—(2.5) 
att = T = 5000 for different values of At = 2. The error is defined as Ey = |v -vy | - |w — uw |, 
where v and w are the numerical solutions for a very fine resolution (At = 0.001), and vy and wy 
are the numerical solutions for larger values of Ат. 


N At EN EN [At 
500 10 0.0923 0.0092 
1000 5 0.0433 0.0087 
5000 1 0.00727 0.0073 
10000 0.5 0.00353 0.0071 
50000 0.1 0.000682 0.0068 


gives At = T/N = 1. In Fig. 2.2, we show the numerical solutions (Un, Wn) in 
a parametric plot, and we note that the solutions are periodic. In electrophysiology, 
such solutions are useful for studying pacemaker cells that keep on creating action 
potentials at a steady rate. 


2.2 What Is the Error? 


In the very simple equation in the previous chapter, we had a formula for the 
analytical solution and a formula for the numerical solution, and thus it was 
straightforward to find the error introduced by replacing a derivative by a difference. 
For the FitzHugh-Nagumo equations, this is harder. In numerical analysis there are 
techniques for proving error bounds for numerical methods. But the proofs tend to 
be very technical and often involve constants that need to be estimated. So we are 
looking for something simpler. If we just assume that we have convergence towards 
the correct solution as Aż tends to zero, we can compute an accurate approximation 
of the solution by using a very small Aż and then monitoring the convergence towards 
this highly resolved solution. 
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Fig. 2.3 Numerical solutions v, (left) and wn (right) of the FitzHugh-Nagumo model specified by 
(2.3)- 2.5). We compare the solutions for a very fine resolution (At = 0.001, solid blue line) to the 
solutions for three cases of coarser resolution (dotted orange line). In the upper panel, we consider 
At = 10, in the middle panel, we consider Ат = 5, and in the lower panel, we consider Ат = 1. To 
improve visibility, we only show the solutions from t = 4000 to t = 5000. As At is decreased, the 
coarse numerical solutions are more similar to the numerical solution computed with a very small 
At, and for At = 1 the two solutions are indistinguishable. 


For simplicity, we assume that we are merely interested in the error at the final 
time = T = 5000. We first compute the solution using an extremely fine resolution 
(At = 0.001) and regard that as the correct solution at time T. Then, we compute 
solutions for varying resolutions (different values of the time step Ar) and compare 
the "correct" and approximate solutions. In Fig. 2.3, we compare the solutions for 
a few different choices of Ar. The error defined by Ey = |v — vy| + |w — wyl, 
where (v, w) denotes the fine scale solution, is given in Table 2.1. Again, we observe 
that Ey /At is more or less constant and we therefore conclude again that the error 
in the numerical solution is proportional to the time step Ar. In other words, the 
convergence is linear (or first order) in At. 
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Fig. 2.4 Illustration of a 
definition of the action 
potential duration, APD50. 
The APDSO value is defined as 
the duration between the two 
times v crosses the threshold 
value defined at the center 
between the maximum and 
minimum value of v. 
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2.3 Upstroke Velocity and Action Potential Duration 


In applications of the FitzHugh-Nagumo model, the unknown function v is often 
used to represent the membrane potential of an excitable cell, e.g., a neuron or a 
cardiac cell, firing a sequence of action potentials. A single action potential from the 
solution of the FitzHugh-Nagumo model is illustrated in Fig. 2.4. In general terms, 
the action potential first consists of a period during which the value of v increases 
slowly, followed by a more rapid increase (the upstroke). Then, v decreases relatively 
slowly for a while before it decreases rapidly back to the minimum value and starts 
increasing again. In this setting, we often refer to v increasing as depolarization, and 
v decreasing as repolarization. We will come back to these terms below where we 
introduce models with proper physical units. 

By solving the FitzHugh-Nagumo model equations for different values of the 
model constants, or parameters, (a, c1, C2, b, and d), we could gain some insight into 
how the parameters affect the firing of action potentials. For example, we could 
investigate how the parameters affect the frequency of firing or the shape of the fired 
action potentials. Two properties that are of interest regarding the shape of the action 
potential are the maximal upstroke velocity and the action potential duration. The 
maximal upstroke velocity is often defined as the maximum value of the derivative 
of v with respect to time. Using a finite difference approximation of the derivative, 
this can be defined as E 
max (==) | (2.10) 


The action potential duration is often defined in terms of a given percentage of 
repolarization, for example APD50 or APD90, for 50% or 90% repolarization, 
respectively. Here, APD50 represents the duration from the start of the action 
potential until the membrane potential reaches a value that is 50% repolarized (i.e., 
at vso = 0.5 (max,(v,) + min,(v,))), at ae Similarly, APD90 is defined as the 
duration from the start of the action potential until the membrane potential reaches a 
value that is 90% repolarized (i.e., at ооо = тах„(о„) — 0.9(max,(v,) — min, (vn ))), 
at poni The start of the action potential used in the definition of the action potential 
duration can, for example, be defined as the point in time when the maximal 
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Fig. 2.5 The numerical solution, v, , of the FitzHugh-Nagumo model, defined by the two equations 
о = civ(v — a)(1 — v) - caw, and ш” = b(v — dw). The parameters are as specified іп (2.5), except 
that in each row of the figure, the value of either a, c1, c2, b, or d is adjusted. The title above each 
plot specifies the parameter change. 


upstroke velocity occurs, or when the membrane potential crosses the 50% or 90% 
repolarization thresholds, 050 or voo, during the upstroke, denoted by ad Or i 
respectively. In the latter case, APD50 and APD90 can be defined as 


APD50 = 199%" — tab, (2.11) 
APD90 = 199%" — рр, (2.12) 


Such a definition of APD50 is illustrated in Fig. 2.4. 

In Fig. 2.5, we show the numerical solution, v, of the FitzHugh-Nagumo model 
with different choices of parameters. In each row, we consider three different values 
of one of the parameters а, су, c», b, or d, and keep the remaining values fixed at 
the values specified in (2.5). In the plots, we observe that increasing the value of 
c, appears to make the action potentials longer and the firing frequency slower, 
whereas the opposite effect is observed when c» or b are increased. In Fig. 2.6, we 
study the effects on the individual action potentials more closely. In the left panel, we 
have zoomed in on the points in time representing the action potential upstroke. We 
observe that decreasing the value of сү reduces the upstroke velocity, but changing 
the other parameters do not seem to have a significant effect on the upstroke. In the 
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Fig. 2.6 The numerical solution, vn, of the FitzHugh-Nagumo model, defined by the two equations 
v' = cıv(v — a)(1 — v) - cow, and ш” = b(v — dw). The parameters are as specified in (2.5), except 
that in each row, the value of either a, сі, с2, b, or d is adjusted. The legends at the right-hand 
side of each row specify the parameter changes. The time axes of the solutions are adjusted such 
that the maximal upstroke velocity occurs at the same time for all the parameter changes. The left 
panel shows the upstroke of the action potential and the right panel shows one action potential for 
each parameter set. 


right panel, we consider a single action potential for the different parameter choices. 
We observe that all the parameters have a significant effect on the action potential 
duration. 
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Chapter 3 wies, 
The Diffusion Equation 


The diffusion equation appears in many applications in science and engineering, 
and computational physiology is no exception. In its most basic form, the diffusion 
equation is also useful as an example of how to deal with a PDE using numerical 
methods. We will start by considering it as a stand-alone model, but in the next 
chapters we will study itin combination with non-linear ODEs. This chapter therefore 
serves as a warm-up for the more complex models. We will also follow the path we 
started above. In the very simplest case of an ODE, we found a formula for the 
solution of both the differential equation and the numerical scheme approximating 
the equation. One nice consequence of this is that, as we saw above, we can explicitly 
study the error introduced by the numerical approximation. In this chapter, we will 
use the same approach to study the error of the numerical approximation of the 
diffusion equation. 


3.1 The Problem = Equation + Initial Values + Boundary 
Conditions 


We consider the diffusion equation 


ди 


д?и 
э”) = 322069. (3.1) 


Неге, и models! the concentration of a chemical in a spatial domain, x is the spatial 
variable and ¢ is time. And since partial derivatives are present, the equation is a 
PDE (see page 13). We study the problem for values of x in the spatial domain (0, 1) 
and for time t in the interval (0, T]. We will assume that the solution is given by u°(x) 
at time f = 0, and that u = О at x = 0 and x = 1 for all time in the interval [0,7]. 


! [n the beginning of these notes, we treat all models as unitless. In models of electrophysiology, 
units can be a true nightmare so we follow the prudent path of trying to deal with one nightmare at 
the time. Units will come later. 
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The function и? = u°(x) is called the initial condition, and u = 0 at x = 0 and x = 1 


are called boundary conditions. In particular, these boundary conditions are referred 
to as Dirichlet boundary conditions, as opposed to Neumann boundary conditions, 
which will be introduced below. 

We will consider the special case of 


иб(х) = sin(zx). (3.2) 

For that particular initial condition, the analytical solution of the diffusion equation 
(3.1) is given by 

u(t,x) = е7"! sin(nx). (3.3) 

This can be verified by observing that the boundary conditions are satisfied since 

clearly u(t,0) = u(t, 1) = 0, since sin(0) = 0, and the initial condition is also satisfied 


since u(0,x) = ѕіп(лх) = u°(x). It remains to be seen whether equation (3.1) is 
satisfied by u(t, x). To this end, we note that 


ди 


m —n2e "^ sin(nx) = —л?и, (3.4) 
and ^ 
д^иб,х) _ Е abu edu (3.5) 
Ox? 


so (3.1) holds. 


3.2 The Numerical Scheme 


We define a numerical approximation of the problem simply by replacing derivatives 
by differences. In order to replace the left-hand side of (3.1), we use the same formula 


as above and state that 
ди u(t + At,x) — u(t, x) 


A 


Ót At 


To approximate the right-hand side of (3.1), we need to recall from calculus that the 
Taylor series of u states that 


(3.6) 


1 
u(t, x + Ах) = u(t, x) + Axux(t,x) + Ах ust, x), (3.7) 


Ou(t,x) д?и(ї,х) А ИК 
where u and uxx are shorthand for ~~ and —777^, respectively. Similarly, 


1 
u(t, x — Ах) = u(t, x) — Axux(t,x) + Ах ихх@, x), (3.8) 


and therefore, by adding (3.7) and (3.8), we get the approximation, 
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u(t, x — Ax) — 2u(t, x) + u(t, x + Ax) 
Ax? | 


Uxx(t,x) © (3.9) 
As above, we introduce t„ = n x At, where At = T/N for a sufficiently large integer 
М. Furthermore, we introduce the spatial mesh points given by x; = (j — 1) X Ax for 
j= 1,..., M, where Ax = 1/(M — 1) fora suitable? integer М. 

We have defined the points ((t,,x;)) and we let и” denote a numerical 
approximation in these points. It remains to define these values. The beginning 
is very simple; we clearly want the numerical scheme to satisfy the initial condition 
so we define, 

u? = u(x;), (3.10) 


and for our special problem we have иў = sin(zxj). Next, we apply the boundary 


conditions and define 
ut = 0, (3.11) 


апа 
uM = 0, (3.12) 


foralln < М. We still need to find values of u7 forj = 2,...,M-landl «n x N. We 
find these by using the finite difference approximations (3.6) and (3.9). Specifically, 
we define и? by the following numerical scheme, 


п+1 п п (0 п п 
uj =u; _ и? 1 2и? и cii 
At Ax? f | 
By defining 
At 

= —, 3.14 
P= 45 (3.14) 

we can write the scheme in computational form, 
we = ри? a + (1 – 2р)и? + QW (3.15) 


Starting with n = 0, we note that we can compute all the values at time f; = At 
since и} only depends on values at the time fo = 0 and they are given by the initial 
condition и? (j = L..., M). When the values at time t; = At have been computed, 
we can use them to compute all the values at t = 2At and so forth. The scheme is 
illustrated in Fig. 3.1. 


3.3 Analytical Solution of the Numerical Scheme and Error 


This subsection is very technical and can be skipped without serious consequences 
(jump to Section 3.4). The purpose of the section is to provide formulas for the 


? "Sufficiently large integer’ and 'suitable integer’ are intentionally vague — stay calm — it will 
become clear (clearer) later on. 
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Unknown solution 
to be computed 


Solution known from 
previous time step 


ім | “ie 


Fig. 3.1 Illustration of the explicit finite difference scheme for the diffusion equation. The solutions 


ue a " и? апа и? -q from the previous time step are used to compute ће new solution, в ,at this 


time step. The scheme is given by (3.15). 


numerical solution and the associated error for a PDE like we did for an ODE in 
(1.13). Such formulas are available only in rare cases, but can be useful for analyzing 
the error of the numerical approximation. 

For a general initial condition, we have to write a program to use the scheme 
(3.15), and we will do that below, but for the special initial condition (3.2), we can 
find a formula for the numerical solution. It is given by 


и? = (1 — Atu)” sinGrxj), (3.16) 


where 
u- 5 sin?(zAx/2). (3.17) 

In order to show that this is the solution, we first check that the formula holds at time 
to. Inserting n = О into (3.16), we get и" = sin(zx;), which matches the numerical 
solution at time £9 = 0 (see (3.2) and (3.10)). Next, we assume that the formula is 
correct at time f,, and show that it also holds at time ¢,,,,, — then the formula holds 
by induction. By inserting (3.16) into (3.15), we get, 
п+1 = 
7 

= (1— Atu)" [psin(zxj1) + (1 2р) sin(ax;)  psin(zxj41)]. (3.19) 


u ри? a +(1 — 2p)uj + ри? (3.18) 


We need two trigonometric identities from calculus to proceed from here, 
1 — cos(y) = 2 sin?(y/2), (3.20) 


and 
sin(x + y) + sin(x — y) = 2cos(y) sin(x). (3.21) 


These identities hold for all values of x and y. We define 
zj = psin(zxj1) + (1 - 2р) ѕіп(лху) + psin(ztxj,1) (3.22) 


and rewrite it like 
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zj = sin(zxj) + p[sin(ztxj-.1) + sintxj41) — 2 sin(zxj)] (3.23) 


sin(zx;) + p[sin(zx; — Ax) + ѕіп(лху + tAx) – 2sin(zx;)], (3.24) 
so by (3.21), we get 


sin(1x;) + p[2 cos(rAx) sin(zx;) — 2 sin(zx;)] (3.25) 
= sin(zx;) + 2p[cos(z Ax) – 1] sin(zx;). (3.26) 


Zj 


By (3.20), we now get, 


zj = sin(ax;) - 4p sin?(zAx/2) sin(zxj) (3.27) 
= ѕіп(лху)[1 — 4p sin? (1Ax/2)], (3.28) 


so, by (3.14) and (3.16), we have 
zj = [1 - Atu]sin(zx;). (3.29) 
Here, и is defined by (3.17). By combining (3.19), (3.22) and (3.29), we find that 
i5! = (1— Arg"! ѕіп(лху), (3.30) 


and thus the formula (3.16) holds by induction. 


3.3.1 What Is the Error? 


Since we now have the analytical solution of the PDE (3.1) given by 
u(t,x) = е7"! sin(zx) (3.31) 
and a formula (3.16) for the numerical solution 


и? = (1 - Ағи)" sin(zx;) (3.32) 


it is straightforward to compare the analytical and numerical solutions. The spatial 
part of the solutions are identical, and therefore it is sufficient to consider the temporal 
part of the solution. We thus consider the relative error in time defined by 


_ | - An)" – | 


e Tta 


En 


(3.33) 


In Table 3.1 we show Ey where N = 1/At for several values of At. We also show 
En /At and we see that, again, this is more or less constant and therefore we again 
have linear convergence with Ey ~ 48 x At. 
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Table3.1 Error ofthe temporal part of the numerical solution of the diffusion equation att = T = 1 
for Ax = 0.001 and different values of At = Z. The error is defined in (3.33). 


N At EN EN [At 
100 0.01 0.406 40.6 
1000 0.001 0.0478 47.8 
10000 0.0001 0.00485 48.5 


100000 0.00001 0.000479 47.9 


3.3.2 More on the Error 


We can go one step further in comparing the analytical (3.31) and numerical (3.32) 
solutions using the formulas for these solutions. Since the spatial dependency of the 
solutions are equal, we just pick x — 1/2 and then we want to compare the analytical 
solution , 

U(T) -uT,1/2 = e* T. (3.34) 


and the numerical solution 
UN = wun = 0 - At)", (3.35) 


at the final time t = T = N x At. Here we have assumed that M is an odd number. 
Recall, again from calculus, that the Taylor series of the sine function states that? 


sin(z) = z + O(z)), (3.36) 


so for small values of z, we have sin(z) « z. We can use this to approximate the value 
of 


m 5 sin (rAx/2) (3.37) 


for small values of Ат. By the Taylor series, we get 
pii eA ud (3.38) 
Ax? ` ` 
So the numerical solution satisfies 
UN x (Ан), (3.39) 


We can also approximate the analytical solution using the following Taylor series for 
the exponential function, 


3 You don’t remember the O-notation? It is shorthand for telling how fast something goes to zero. 
So f(x) = О(х?) means that f(x) goes to zero as fast as x? goes to zero for small values of x. In 
numerical analysis this notation is often used to indicate the size of an error term. Often, it is the 
remainder of a truncated Taylor series. 
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е4 =1-:+ 0(22). (3.40) 


By using this approximation, we find that the analytical solution at x = 1/2 сап be 
approximated by 


U(T) = e 7T = e NN = (g7 AN с (р AD, (3.41) 


hence, clearly, we have 
UN x U(T). (3.42) 


3.4 Instabilities in the Numerical Solution 


Numerical instabilities often appear and it can be hard to understand the origin. Here, 
we will show instabilities that appear because of too long time steps, but keep in 
mind that a long list of other instabilities can appear as well. Differential equations 
can have oscillatory solutions, but as the mesh parameters are refined, the numerical 
solution should converge towards the correct solution. Numerical instabilities, on 
the other hand, manifest themselves by diverging as the mesh is refined, rather than 
converging. 


3.4.1 Specification of a Stable Problem with Unstable Numerical 
Solutions 


Suppose you have a tank of length 1 filled with water, where you have ink added 
to the water for x x 1/2 but no ink for x > 1/2. In the middle of the tank, at 
x = 1/2, there is an impermeable membrane, so there is no ink leaking from the left 
to the right side of the tank. At time т = 0, we remove the membrane and we are 
curious about what is going to happen. Intuitively, we expect the ink to diffuse to the 
right-hand side of the tank and that, eventually, the ink will be uniformly distributed 
throughout the tank. Qualitatively, we can get an impression of what happens by 
solving the diffusion equation 


ди д?и 
ap) = gato (3.43) 


We will use the initial condition u(0, x) = 1 for x < 1/2 and u(0, x) = 0 for x > 1/2. 
The boundary conditions are given 


ðu ðu 
5, 60) = ax 1) =0, (3.44) 
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i.e., that the spatial derivative is zero at the boundary. This is referred to as a no-flux 
boundary condition, or a Neumann boundary condition, and simply means that we 
don't allow the ink to leak out of the tank. 


3.4.2 Numerical Scheme for Neumann Boundary Conditions 


In order to deal with this alternative set of boundary conditions in the numerical 
scheme, we consider the two Taylor series of u considered in Section 3.2, i.e., 


1 
u(t, x + Ax) = u(t, x) + Axu,(t,x) + 5x utt, x), (3.45) 


and 
1 
u(t, x — Ax) = u(t, x) — Axu,(t,x) + Ах ux). (3.46) 
By subtracting (3.46) from (3.45), we obtain 
u(t, x + Ax) — u(t,x — Ax) = 2Axus(t, x), (3.47) 


which yields 
u(t, x + Ax) — u(t, x — Ax) 


2Ax 
Replacing the derivatives by this difference in the boundary condition (4.24), we get 


(3.48) 


Ux(t,x) = 


uy — Uo Wy ИМ-1 
= a = E 
2Ax 0, 2Ax E pun 
which yield 
ш = U3, им+ = Им-1› (3.50) 


for all n x М. Inserting (3.50) into the main scheme (3.15) for j = 1 and j = M, we 
obtain 


utt! = (1-2p)u! + 2pu5, (3.51) 
e -(1- 2p)u, + 2риђ 1: (3.52) 


3.4.3 Example of Instabilities in the Numerical Solution 


In Fig. 3.2 we have solved the problem numerically using At = 0.0001 and Ax = 0.02 
and we note that the solution seems to evolve as we expected. But in Fig. 3.3, we 
attempt to take longer time steps, and then we see pretty wild oscillations in the 
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numerical solution. This is a classical type of numerical instability that often appears 
in numerical methods. Usually it helps to reduce the time steps, and it certainly helps 
in this case. 


ї=0 t - 0.01 t2 0.05 t- 0.1 = 0.5 
1 1 1 1 1 
3 05 05 оз ааз 
0 0 0 0 0 
0 05 10 05 1! 0 05 10 05 1 0 05 1 
х х X X X 


Fig. 3.2 Numerical solution of the diffusion equation with boundary conditions ди (1, 0) = 
ди (т, 1) = 0 and initial conditions и(0, x) = 1 for x < 1/2 and u(0, x) = 0 for x > 1/2 
using Ат = 0.0001 and Ax = 0.02. The solution is plotted at five different time points. 


t-0 
1 
5 0.5 
0 
0 0.2 0.4 0.6 0.8 1 
t = 0.001 
2 
> 0 
-1 
0 0.2 04 0.6 0.8 1 
x 10? = 0.01 
Fig. 3.3 Numerical solution A 
a | 


of the diffusion equation 


with boundary conditions -2 

ди (т, 0) = ди (1, 1) = 0 and 0 0.2 0.4 0.6 0.8 1 
initial conditions u(0, х) = 1 2X 10% t= 0.05 

for x € 1/2 and u(0, x) = 0 

for x » 1/2 using At — 0.001 й: 

and Ах = 0.02. For these Es 

values of At and Ax, we get 0 0.2 0.4 0.6 0.8 1 
pretty wild oscillations in 1 x 109^ 1=01 

the numerical solution. Note 

that the scaling of the y-axis 3.0 —À Ne 
is different in each row of -1 

the figure, corresponding to 0 02 04 06 08 1 


different points in time. x 
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3.5 Stability Condition for the Numerical Scheme 


The diffusion equation satisfies a maximum principle stating that the solution will 
always be bounded by the values of the initial condition and the boundary conditions. 
We will show that the numerical solution generated by the scheme (3.15) satisfies 
the same principle, provided that the time step is sufficiently small. For simplicity 
we again assume that the boundary conditions are given by и = 0 for x = 0 and 
x = 1. Concretely, we consider the scheme 


ил! = ри? 1 + (1 – 2р)и? + ри?» (3.53) 
where we recall that А 
t 
= —, 3.54 
A (3.54) 
and where the boundary conditions are given by и] = 0 and иу, = 0. Define 
u, = max |u®], (3.55) 
j 


i.e., the biggest (in absolute value) initial value, and assume that p < 1/2, i.e. we 
assume that 


А 2 
Аа, (3.56) 
2 
Now, we want to prove that 
ГА < u4 (3.57) 


forall j € [2, M - 1] and n > 0. We will show this by induction and start by assuming 

that (3.57) holds for an arbitrary value of n. Then, by the scheme (3.53), we get* 
| = [рит , + (1 7 2р)и? + ри", | 

< ош" |+ Q1 - 2p)lu?] + plu, 

< pus + (1 — 2p)u. + pus 


= Uy 


and therefore (3.57) holds by induction. 

We also note that the criterion (3.56) is in agreement with the results observed in 
the numerical example in the previous section. In that case, we had m = 0.0002. 
Thus, for the stable case of Ат = 0.0001 (Fig. 3.2), the stability criterion (3.56) 
was satisfied, whereas for the unstable case of At — 0.001 (Fig. 3.3), the stability 
criterion was not satisfied. 

When using explicit finite difference schemes to solve equations where the 
diffusion equation is part of the problem, a stability condition similar to (3.56) 
usually has to be satisfied in order to obtain proper numerical results. 


^ Here we use the triangle equality stating that |a + b| < |а| + |b] for any numbers a and b. 
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3.6 A Brief Comment on Uniqueness 


We have seen that a formula can be obtained for the solution of the diffusion equation 
and for the numerical approximation of the same problem. But are these solutions 
unique? Can there be other solutions than those given by the formulas? In order to 
see that the solution of the continuous problem is in fact unique, we assume that we 
have two solutions, и and v, with coinciding initial condition given by f = f(x), and 
the usual Dirichlet boundary conditions. Since u and v solve 


ди д?и 
a, 09 = zat” (3.58) 


and 3 
Ov д^ 
5; 09 = gato (3.59) 


respectively, we find (by subtracting (3.59) from (3.58)) that the difference between 
the solutions, given by e = u — v, solves the equation 


де д?е 
ar x)= ag x). (3.60) 


Now, we can define a measure of the difference between these solutions as a function 
of time, 


ШЗ! 
E(t) = = i e7(t,x)dx, (3.61) 
2 Jo 
and observe that 
1 1 2 
E'(t) = | ie” dos f e(t, 1E). (3.62) 
0 ді 0 дх 
By using integrations by parts, we find that 
! д?е ! (8e : 
/ ex) oen) = -f (Ze) dx, (3.63) 
and therefore, 
E'(t) x 0. (3.64) 


Since и and v have the same initial condition given by f = f(x), we clearly have 
е(0, x) = 0 for all relevant values of x. Therefore, since E(0) = 0, and E’(t) x 0 
we have E(t) = О for all time and thus и and v are equal, and the solution of the 
problem must be unique. A similar argument can be given for the discrete case, so 
the numerical solution given by the formula (3.16) is also unique. This is a classical 
energy argument and it can be extended to also provide stability estimates of the 
solutions; see, e.g., [1]. 
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A 
Check for 


Chapter 4 шашы 
Implicit Numerical Methods 


In the previous chapter, we saw that the simple explicit numerical scheme resulted 
in an instability problem. We also saw that the problem could be resolved by 
using sufficiently short time steps. But in many situations, short time steps become 
exceedingly short, as can be seen, e.g., in the stability criterion (3.56). This means 
that we have to perform computations for a very large number of time steps to reach 
the final time and it is therefore tempting to look for alternatives. The most common 
alternative is to use an implicit scheme, which generally allows for much longer time 
steps. 


4.1 Explicit and Implicit Numerical Schemes 


In Section 1.4 (page 7) we briefly introduced the concept of explicit and implicit 
numerical schemes. So far we have only considered explicit schemes and the reason 
for that is just simplicity. If we have a differential equation that can be written in the 
form 

u'(t) = F(u(t)), (4.1) 


we can, as seen in Chapter 1, use the the approximation 


u(t + At) — u(t) 


(t) = 4.2 
u(t) x (4.2) 
to replace (4.1) by the difference equation, 
yor! = у" 
—— = F(u”), 4.3 
A = Еи") (4.3) 
leading to the explicit scheme 
ult! = и" + AtF(u"). (4.4) 
© The Author(s) 2023 33 
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But an alternative, equally plausible, approach is to replace (4.1) by the difference 
equation 


yor! — и" 
—— = F(u™!), 4.5 
Ap (и!) (4.5) 
which leads to the implicit scheme 
u”! — AtF(u”*!) = u”. (4.6) 


As we shall see below, the advantage of the implicit scheme is that we obtain 
numerically stable results for longer time steps, but the disadvantage is that we have 
to solve a potentially nonlinear equation of the form 


u— AtF (u) = й (4.7) 


at every time step. Here, iiis known from the previous time step, and u is the unknown 
that we need to compute. 

We will now show how to derive implicit schemes and demonstrate that they are 
more stable than the explicit versions. But keep in mind that we also have to deal 
with accuracy. For this purpose we still want the time steps to be reasonably short 
without being too short. 


4.2 An Implicit Scheme for the Diffusion Equation 


Let us first recall how we derived the explicit scheme for the diffusion equation!, 
Ut = Uyy, (4.8) 
equipped with the initial and boundary conditions, 


u(0, x) = и°(х), (4.9) 
u(t,0) = u(t, 1) = 0. (4.10) 


As in (4.3), we can approximate this equation by replacing the time derivative by 


u(t + At) — u(t) 
At ` 


Ut © (4.11) 
For the right-hand side of (4.8) there are two alternatives (actually many alternatives). 
We can either approximate the right-hand side by a difference approximation of 
Uxx(t, x) or of uxx(t + At, x). The first alternative results in the explicit scheme we 
used above, 


! [f you just browse these notes and didn't recognize the notation used here, we repeat that 
и = Ou/Ot, and uxx = uJàx? — keep on browsing. 
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n+l _ „п n _ n n 
и? иу Ид 2и? +? a2 
А Ах? , $99) 
and the second alternative results in the implicit scheme, 
n+l уп n+l . п+1 n+l 
u; иу Mj 2и? + ит NT 
At Ах? | SM 


In Fig. 4.1 we have illustrated how this scheme works. Here, we realize that there 
is a fundamental difference between the explicit (see Fig. 3.1) and implicit schemes. 
The explicit schemes are very simple since every value is simply an explicit function 
of the values of the previous time step. This is straightforward to implement in 
software. But the implicit scheme is much more convoluted. In fact, all values at 
time f,4; are coupled with all other other values at that the same time step. Well, 
actually, every point is connected to two neighbors, but these again are coupled to 
their neighbors and so on. So we end up with a linear system of equations. 


Unknown solution 
to be computed 


б з а == Solution known from 
n previous time step 


Fig.4.1 Illustration of the implicit finite difference scheme for the diffusion equation. The unknown 
solutions ur andu ai from the present time step, in addition to the known solution и? from the 
previous time step, are all needed to compute the solution ut. The scheme is given by (4.14), or 
on matrix form by (4.21). 
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We will write the scheme (4.13) as a linear system of equations, but first we rearrange 
it by putting unknowns (variables at time tn+1) at the left-hand side and previously 
computed variables (ї„) on the right-hand side of the equation, 

n+l 


- pui, + (1+ 2p"! — рит =u}, (4.14) 


where again 


At 
p= 5: (4.15) 
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The implicit scheme has the form (4.14) for j = 3,..., M — 2 but takes a different 
form for j = 2 and j = M - 1. Recall that the boundary conditions are specified by 
(4.10), given by 

uy = иу = 0 (4.16) 


for all n. Hence, for j = 2, ће scheme (4.14) takes the form, 
(1+ 2p)u5*! — рит! = ит. (4.17) 


Similarly, for j = М — 1, we get the scheme 


- ри!) + (1 +2р)и]! = ит. (4.18) 


We are now ready to rewrite the somewhat messy scheme defined by (4.14), (4.17), 
(4.18) in matrix form. To this end, we defined the vectors? 


u, 
us 
u"=| ; |, (4.19) 
VAM 
and the matrix 
1+2р -p 0 € 0 
-p 1+2р -p 0... 0 
А=| 0 p 1+2p-p 0: 0 |. (4.20) 
0 ee 0 -p 1+2p 


With these definitions, we can rewrite the scheme defined by (4.14), (4.17), (4.18) 
on the simple form 

Au”! = и". (4.21) 
The matrix A defined in (4.20) is tridiagonal and this makes the system (4.21) easy 
to solve. In the software associated these notes, the solution is shown in Matlab, but 
the system is easy to solve in any numerically oriented computing system. 


4.3.1 Neumann Boundary Conditions 


In Section 3.4 we considered the numerical solution of a specific diffusion equation 
problem given by 


? Note that uï and и», are given directly by the boundary conditions at every time step, so these 
values do not have to be included in the vector of unknowns. 
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Ut = yx. (4.22) 


with the initial and boundary conditions 


1, 1Ёх< 1/2, 
u(0, x) = | (4.23) 
0, їЁх> 1/2, 
ux(t,0) = ux(t, 1) = 0. (4.24) 
For the inner points j = 2,..., М — 1, an implicit numerical scheme for this problem 


can be given by (4.14) like above, but for the boundary points (j = 1 and j = M), we 
need to take the different set of boundary conditions into account. One way to do this 
is by following the same procedure as in Section 3.4.2. That is, using the difference 


u(t, x + Ax) — u(t, x — Ax) 


SA (4.25) 


ux(t, x) = 


Replacing the derivatives by this difference in the boundary condition (4.24), we get 


uy — цу има 7 Ug 
= 0), ————7— =), 4.26 
2Ах 2Ах ( ) 
which yield 
ш = U5, UM+1 = Им-1› (4.27) 


for all n x М. Inserting (4.27) into the main scheme (4.14) for j = 1 and j = M, we 
obtain 


(1+ 2put*! – 2pu5*! = и", (4.28) 
(1+ 2p)uM! – 2рип+! = и", (4.29) 
and the scheme can be written оп matrix form 
Au"! = и", (4.30) 
where 
uy 
u” 
и =| 2 |, (4.31) 
и" 
апа 
1+2p -2p 0 0 
-p 1+2p -p 0 0 


А=| 0 -p 1+2р-рО.. 0 J (4.32) 
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4.4 The Implicit Scheme Is Stable 


We noticed above that if the time steps we used for the explicit scheme were too long, 
we got a solution with wild oscillations. In Fig. 4.2 we repeat these computations 
using the implicit scheme (4.30)-(4.32) and we notice that the oscillations are gone. 
The computations could indicate that the solutions are stable for any value of At. 
In fact, unconditional stability of the implicit scheme for the diffusion equation is 
classical and can be proved. One argument can be found in [3]. 


Explicit scheme А Implicit scheme 


x 1046 


t2 0.1 t=0.1 
> 0 505 mE a 
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Fig. 4.2 Numerical solutions of the diffusion equation with boundary conditions ди (1, 0) = 
ди (т, 1) = 0 and initial conditions u(0, x) = 1 for x < 1/2 and и(0, x) = 0 for x > 1/2 
using Ат = 0.001 and Ax = 0.02. The left panel shows the solution for the explicit scheme (as also 
seen in Fig. 3.3), and we observe wild oscillations. In the right panel, we show the solution of the 
implicit scheme, and the solution appears to be more reasonable. Note that at T = 0.5, the explicit 
scheme is so broken down that the solution returned by the scheme are not numbers (NaNs). 
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4.5 Comments and Further Reading 


1. As mentioned above, when the PDE under consideration is in one spatial 
dimension, the associated linear system (4.21) is easy to solve. But if we consider 
the diffusion equation in three spatial dimensions, e.g., 


Ut = Uxx + Uyy + Uzz, (4.33) 


this problem becomes much harder. The solution of linear systems in the form of 
Ax = b, where A is a matrix, b is a known vector, and x is unknown, is a crucial 
problem in scientific computing. It is a fundamental part of almost all scientific 
computing projects and presents significant challenges. While small systems are 
relatively straightforward to solve, larger systems become increasingly difficult. 
This problem is well studied in the field of scientific computing, but it remains 
a challenging task in general. 

2. Numerical linear algebra is a field of research concerned with operations on 
vectors and matrices on computers. A thorough introduction to the field is given 
in [1]. 

3. A classical introduction to iterative methods for solving linear systems is 
presented in [2]. 
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A 
Check for 


Chapter 5 peas 
Improved Accuracy 


In the examples we have considered so far, the error introduced by the numerical 
scheme has been O(At). This can easily be improved as we will show through two 
simple examples. 


5.1 Back to the Simplest Equation 


The first equation we studied was 


y =y, (5.1) 
and we approximated the equation using the finite difference approximation 


Yn+1 = Yn 


=з; 2 
N y (5.2) 


This scheme resulted in an error proportional with the time step At. From Chapter 4, 
we realize that we can also use the approximation 


Yn+1 = Yn 


до“? (5.3) 


And, as a compromise between these two alternatives, we can use the following 
midpoint! approximation, 


Уп+1 = Yn 


1 
At = Оп + Уп+1). (5.4) 


The explicit, implicit and midpoint schemes can be written on computational form 
as follows, 


! This scheme is sometimes called the midpoint scheme (for obvious reasons) and sometimes called 
the Cranck-Nicolson scheme because it was developed by John Crank and Phyllis Nicolson, [1]. 
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Уп+1 = (1 2n А)у», (5.5) 
1 
Уп+1 = 1- AP" (5.6) 
1+ At/2 
ntl = тт Уп 7 
ун = ayy) (57) 


respectively. In Table 5.1 we show the errors introduced by these three schemes 
and we note that the midpoint scheme is clearly more accurate than the two other 
schemes. Furthermore, we find that the error of the implicit and explicit schemes 
are both proportional to At, or O(At), whereas the error of the midpoint scheme is 
proportional to Ax?, or O(At?). This means that the implicit and explicit schemes 
have first order (linear) convergence, whereas the midpoint scheme has second order 
(quadratic) convergence. 


Table 5.1 Errors of numerical solutions of the differential equation (5.1) with initial condition 
uo = latt = Т = 1 for different values of Ar. The errors are defined as Ez = |y(1) – ум, |. 
Е; = |y(1) - yu i| and Em = |y(1) – ум, т |, where y(1) = e is the analytical solution, and yy e, 
Ум i; and ум, т, are the numerical solutions of the explicit (5.5), implicit (5.6) and midpoint (5.7) 
schemes, respectively, at time 1. 


At Ee Е„ /Аї Е; Ei | Nt Em Е,„ /At? 
0.1 0.125 1.25 0.15 1.5 0.00227 0.227 
0.01 0.0135 1.35 0.0137 1.37 2271-10? 0.227 
0.001 0.00136 1.36 0.00136 1.36 227-1077 0.227 
0.0001 0.000136 1.36 0.000136 1.36 227-107? 0.227 


The difference in numerical errors reported in Table 5.1 may not seem to be 
dramatic. But if you require the error to be less than, say, 10-10, then the number of 
time steps needed for the first order schemes are about 2.85 x 10° larger than what is 
needed for the second order scheme. So, the difference in computing time between a 
first and second order scheme can be dramatic. For the very simple model considered 
here, the computation is trivial in any case, but with a challenging system of partial 
differential equations in three spatial dimensions, the difference in computing efforts 
can become enormous. 


5.2 A Linear Reaction-Diffusion Equation 


We considered the diffusion equation above. That equation becomes a little more 
interesting if we add a reaction term to it, 


Uy = Uxx + f(u). (5.8) 
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Here, f — f(u) is referred to as a reaction term. In order to keep things simple, we 
will consider a linear reaction term and define 


f(u) = au, (5.9) 


and we will wait a little before defining the value of Д. Above, we noted that a 
numerical scheme for the diffusion equation could be written in a very compact form 
by introducing vector/matrix notation. We will extend this in order to define three 
alternative schemes for the reaction-diffusion equation (5.9). Note that we still apply 
the following initial and boundary conditions, 


и(0, x) = u(x), (5.10) 
и(1,0) = u(t, 1) = 0. (5.11) 
As above, we let 
ил 
из 
u” = |, (5.12) 
Wy a 


and in addition, we introduce the matrix 


—-2] (о 0 
АКЕЛУ С. 

р=-—; И RE (5.13) 
(ЖЕ 01-2 


With this notation at hand we can approximate the linear version of equation (5.8) 
using an explicit, implicit or midpoint approximation of the right-hand side of the 
equation 


n+l _ „п 

: x “= Du" + Àu", (5.14) 
n+l n 

E x 2 = Du! + aunt, (5.15) 
n+l _ „п 1 

= N е ;0Dw + Du! + Au” + aw), (5.16) 


To help ease the implementation of the numerical schemes and make it more 
straightforward to write the schemes in matrix form as in, e.g., (4.21), it is convenient 
to write these schemes in computational form, that is, collect all terms involving the 
unknown solutions to be computed in each time step, u"*!. on the left-hand side and 
all the terms involving the know solutions from the previous time step, и”, on the 
right-hand side. The schemes (5.14)-(5.16) can be written in this form by, 
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= [(1 + AADI  AcD]u" (5.17) 

(1 -AAt – AtD]u™! = u” (5.18) 
AAt AAt At 

ТЕТ? (+) в 


where / is the identity? matrix. 
Let us now consider the problem (5.8)-(5.11) with a specific choice of A = 2л? 
and u°(x) = ѕіп(лх). Then, the analytical solution is given by 


u(t,x) = е"! sin(zx). (5.20) 


In Fig. 5.1, we show the results of the implicit (5.18) and midpoint (5.19) schemes 
together with the analytical solution at time T = 1. In addition, we show the results 
of the explicit scheme at T = 0.1. We have used Ax = 0.01 and Ar = 0.001. For this 
choice of discretization parameters, we get wild oscillations for the explicit scheme. 
The implicit and midpoint schemes produce more reasonable solutions, and we note 
that the midpoint scheme is clearly more accurate than the implicit scheme. 
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Fig. 5.1 Analytical and numerical solutions of the reaction-diffusion equation (5.8) with A = 2л?, 
boundary conditions u(t, 0) = u(t, 1) = О and initial conditions u(0, x) = ѕіп(лх). In the numerical 
schemes, we use At = 0.001 and Ax = 0.01. In the upper panel, we show the solution of the implicit 
and midpoint schemes together with the analytical solution at time T = 1. The lower panel shows 
the solution for the explicit scheme at T = 0.1, and we observe wild oscillations. 
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A 
Check for 


Chapter 6 C 
A Simple Cable Equation 


The cable equation was first derived to model transport of electrical signals in 
telegraphic cables. But it later gained enormous popularity as a model of transport 
of electrical signals along a neuronal axon. In Chapter 9, we will discuss how this 
equation is derived and how the different terms in the equation come about. But 
here, we will just take a simple version of the equations for granted and then try to 
solve them. We will observe that the few techniques we learned above are actually 
sufficient to solve the non-linear reaction-diffusion equations we consider here. 


6.1 A Non-Linear Reaction-Diffusion System 


We will consider a system of equations where we add a diffusion term to the 
FitzHugh-Nagumo equations. Specifically, we consider the equations!, 


ш = охх + с10(о — a)(1 — v) — cow, (6.1) 
ш, = b(v — dw). (6.2) 


We use the same constants as above, 
a = —0.12, cy = 0.175, c; = 0.03, b = 0.011, d = 0.55, (6.3) 


and in addition we use the diffusion coefficient? 6 = 5-107. 


! Recall, once again, that we use the convention that v; = бе and охх = дь : 
2 We still use no units; they will be introduced later. 
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6.2 The Explicit Numerical Scheme 


We note that we can use the same approximation of о, and ш; that we used for 
the FitzHugh-Nagumo model in Chapter 2 (see (2.6) and (2.7)), and the same 
approximation for v,, as we used for the diffusion equation in Chapter 3 (see (3.13)) 
to define an explicit numerical scheme, 


vr = и" о" 1 — 20" + Uf 
ee F187 On - 401 — vj) - cw), (6.4) 
nel _ ш" 
— Í 2 p(y" — аш") (6.5) 
At 7 JU | 


Here, v^ and ш” denote approximations of v(x;,tn) and w(x;,tn), respectively. It is 
straightforward to put this scheme in computational form, 


п+1 
" 


п+1 _ п n n 
ш =W; * Atb(v; dw;) 


= ро 1 + (1 2р)о? + Роў + At[cio7 (Un - ay1- vj) - с2ш? |, 


where р = бАт/Ах?. 


6.3 Traveling Wave Solutions 


Traveling wave solutions are characteristic of solutions of reaction-diffusion models 
of neuronal axons and myocardial tissue. Here, we will see such solutions for the 
simple model given by (6.1) and (6.2). 

In Fig. 6.1, we show the numerical solution of v as a function of x at five different 
points in time. In order to initiate a wave traveling from left to right, we let the initial 
conditions be specified by оо = 0 and wo = 0 for all values of x, except that we set 
vo = 0.26 for x < 0.04. We use Ax = 0.01 and At = 0.005. The boundary conditions 
for v are given by 2 (1,0) = Hr, 1) 2 0. 

In the leftmost panel of Fig. 6.1, at t£ = 0, we see the initial conditions for the 
variable v. In the second panel, at t = 50, we see that the value of v has started to 
increase in an area close to the left boundary of the domain, and in the next panels, 
at t = 100, ғ = 200, and т = 300, we see that the increase in v gradually moves from 
left to right like a traveling wave. 


6.3.1 Adjusting Parameters 


In Fig. 6.2, we show the solution of the model at five different points in time for 
three different values of the parameters ó. We observe that for a low value of ó 
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Fig. 6.1 Numerical solution of v in the FitzHugh-Nagumo equations with a diffusion term added. 
The boundary conditions are given by до (т, 0) = де (t, 1) = 0, and the initial conditions are 


w(0, x) = 0 everywhere and v(0, x) = 0 for x > 0.04 and 0(0, x) = 0.26 for x < 0.04. We use 
At = 0.005 and Ax = 0.01, and show the solution at five different points in time. 


(5 = 1- 1075), the traveling wave moves more slowly through the domain. For 
example, at т = 200, the wave has crossed about half of the distance from x = 0 to 
x = 1 for the default case of 6 = 5 · 1075, and only about a quarter of the distance for 
ó = 1- 1077. Moreover, for a high value of 6 (ô = 20 1075), the wave moves more 
quickly through the domain, and has almost crossed the entire domain at t = 200. 
We also observe that the slope of the wavefront of the traveling wave appears to 
become less steep at the value of 6 is increased. 

In Fig. 6.3, we similarly consider the traveling wave solutions for three different 
values of the parameter cj. As when we adjust the value of б, we observe that the 
wave moves more quickly as the value of c, is increased. 


6.3.2 Conduction Velocity 


In Fig. 6.2 and Fig. 6.3, we observed that the traveling wave moves more quickly as we 
increased the value of б or сі. The conduction velocity is often used to characterize 
the speed with which a traveling wave traverses the domain. In Fig. 6.3, we have 
computed the conduction velocity for some different values of б and cı. Here, we 
have defined the conduction velocity as 


. Ada edu 


CV з 
h-t 


(6.6) 
where x, = 0.5 and x2 = 0.7. Furthermore, tı and t are the points in time when 
the value of v first increases to a value v > 0.5, in the spatial points x; and x9, 
respectively. As expected Fig. 6.3 shows that increasing 6 or c in the model increases 
the computed conduction velocity. 
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Fig. 6.2 Numerical solution of о in the FitzHugh-Nagumo equations with an added diffusion 


term at for three different values of б. The remaining parameter values are as specified in (6.3). 
The boundary conditions are given by 22(r, 0) = 22 (t, 1) = 0, and the initial conditions are 


w(0, x) = 0 everywhere and 0(0, x) = 0 for x > 0.04 and v(0, x) = 0.26 for x < 0.04. We use 
At = 0.005 and Ax = 0.01, and show the solution at five different points in time. 
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Fig. 6.3 Numerical solution of v in the FitzHugh-Nagumo equations with an added diffusion term 
for three different values of 6 for three different values of cı. The remaining parameter values are as 
specified in (6.3) and 6 = 5 · 1075. The boundary conditions are given by i, 0) = др (1, 1) = 0, 
and the initial conditions аге ш(0, x) = 0 everywhere and 0(0, x) = 0 for x > 0.04 and 0(0, x) = 
0.26 for x < 0.04. We use At = 0.005 and Ax = 0.01, and show the solution at five different points 
in time. 
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Fig. 6.4 Conduction velocity for different values of 6 and c1, computed from numerical solutions 
of the FitzHugh-Nagumo equations with an added diffusion equation term. The conduction velocity 
is computed as defined in (6.6). 
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Chapter 7 S" 
Operator Splitting 


In mathematics, a common approach to solving a new problem is to break down 
the problem into sub-problems that you know how to solve and then try to glue the 
pieces together to yield a solution of the new problem. This approach is also very 
useful in software development; well-tested pieces of software can be glued together 
in order to obtain solutions to a wider class of problems. Operator splitting is a 
technique that illustrates this principle very well. Rather complex equations can be 
broken down into more familiar problems and solved individually. It's a miracle that 
it works, but it does. And it's no miracle because there are proofs of convergence. 
Anyway, we will illustrate operator splitting with two examples and then come back 
to this technique when the equations become more challenging. 


7.1 Numerical Schemes for a Reaction-Diffusion Equation 


Suppose we want to solve the following reaction-diffusion equation, 
Ut = Uxx + f(u), (7.1) 
with initial and boundary conditions, 


u(0, x) = u°(x), (7.2) 
u(t,0) = u(t, 1) = 0. (7.3) 


Clearly, we can use the techniques introduced above. By using the notation introduced 
on page 43, we can write explicit, implicit, and midpoint schemes as follows, 
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n+l _ „п 


и и n n 
ae = Du" + flu"), (7.4) 
п+і „п 
u © u = Ри"! +f (a), (7.5) 
yet 1 п п+1 1 п п+1 
xL = Xu + Dut) 2 (f?) + Fu"), (7.6) 


where we have tacitly extended f to be a vector valued function with components 
fi = f (ui). These schemes can be rearranged to computational form as follows, 


и" = [I + A(D]u? + At f(u”), 
= AtD]u”*! — Arf) = u”, 


At п+1 At п+1 At n At n 
\ D 5 Flu = (r+ Sp) un ros 
where, again, we use the convention that known quantities are on the right-hand side 
of the equations and the unknowns are at the left-hand side. The explicit scheme is 
straightforward to implement since computing u”*! simply amounts to evaluating 
the right-hand side. But the implicit scheme has become more complicated than we 
are used to because we have to solve a potentially non-linear system of algebraic 
equations in order to compute и"+!. We can do that using Newton’s method, but we 
can also avoid this by introduction operator splitting. 


7.2 Operator Splitting for a Reaction-Diffusion Equation 


Let u” denote an approximation of u(t,,x), where tn = nAt as usual. Then, the 
problem (7.1) can be solved by alternately solving и; = uxx and и; = f(u). More 
precisely, we assume that an approximate solution has been computed for time f = ty. 
Then, an approximate solution of (7.1) at їл+1 can be found in two steps. First we 
solve 

Ит = Uxx, (7.7) 


from f, to fn+1 with the initial condition u(t,,-) = u” and boundary conditions given 
by (7.2) and (7.3). We let u"*!/2 denote the solution of the first step. In the second 
step, we solve 

и, = f(u), (7.8) 


using the initial condition u(t,,-) = u"*'?. Now, we have broken the somewhat 


complex problem (7.1) down to two problems we are more familiar with. We will 
show how this works with a couple of examples. Note, however, that the main 
message from this chapter is the technique of breaking down the numerical solution 
of a problem into two simpler problems like described in this subsection. So if you 
reach a point where you feel that the mathematics in the remaining part of this chapter 
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Fig. 7.1 Numerical solution 
of the differential equation 6 
и, = Uxx — u’, with boundary 
conditions u(t, 0) = u(t, 1) = 
O and initial condition 
u(0,x) = sin(zx). The 
problem is solved using the 
operator splitting procedure 
described in (7.10)-(7.13) 
with Ax = 0.01 and 

At = 0.001. 


x10? 


becomes a bit overwhelming, it might be good to jump ahead to Chapter 8, where 
we will start applying the methods introduced above to models of electrophysiology. 


7.2.1 Numerical Example 


In the first numerical example, we consider the problem 
_ 2 
Ut = Uxx ти, (7.9) 


with the boundary conditions u(t,0) = u(t, 1) = 0 and the initial condition u(0, x) = 
sin(zx). We want to use the operator splitting procedure introduced above to solve 
this problem numerically, and this will result in two steps. The first step is to solve 


Ит = Uxx, (7.10) 
with an implicit scheme (see page 34). The second step is to solve 
и; = –и?, (7.11) 


using an implicit scheme. Note that this equation has to be solved for each mesh 
point х;, so the implicit scheme reads 


41/2 
unt c ntt 


== == = (uy, (7.12) 


where и" 1/2 is the result of the first step. By using operator splitting we have avoided 
solving a big system of non-linear equations. Instead we need to solve a usual linear 
system of equations arising from the discrete version of (7.10) and a series of second 
order algebraic equations given by (7.12) whose solutions are given by the quadratic 
formula, 
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eda eae 
п+1 J 


- | 7 
“j Y (7.13) 


In Fig. 7.1 we show the solution of this problem using Ax = 0.01 and Aż = 0.001. 
Furthermore, in Table 7.1 we have computed the error by comparing the solution to 
a very fine scale solution of this problem using a standard explicit scheme. Again we 
note that the error seems to be first order in At. 


Table 7.1 Errors of the numerical solutions of the differential equation u; = uxx — u2, with 
boundary conditions u(t,0) = u(t,1) = 0 and initial condition u(0, x) = sin(ztx) for different 
values of At. The error is defined as E = max; |ue,; — uN |, where ue, ; is the solution of the 
problem found using a standard explicit scheme with a very fine time step (At = 1079), andu is ће 
numerical solution of the operator splitting scheme described in (7.10)-(7.13). We use Ax = 0.01. 


At E E/At 

0.01 2.77-1075 0.0028 
0.005 1.27-1075 0.0025 
0.001 2.37-10—% 0.0024 
0.0005 1.17-1076 0.0023 
0.0001 2.35.1077 0.0024 


7.2.2 А Detour via Numerical Integration 


Very few things come for free, but second order convergence does. With minimal 
change of the algorithm above, we can obtain second order convergence. Actually, 
just the first and the last step of the algorithm are slightly changed (half step instead 
of whole step). Exactly the same alteration is present in changing from a first order 
to a second order scheme of numerical integration. We will show that similarity here 
because it might help make the second order operator splitting algorithm seem less 
mysterious, but if you are not interested, just skip this subsection and go to page 58 
where you can read about second order operator splitting. 


First and Second Order Schemes for Numerical Integration 


Suppose we have a smooth function g = g(x) defined on the interval from 0 to 1, 
and that we want to compute an approximation of the integral of g on this interval. 
We start by defining Ax = 1/(M — 1) and х; = (i — 1) x Ax for a sufficiently large 
integer M. In Fig. 7.2, we show two approximations of the function g. The first 
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Piecewise constant Piecewise linear 
approximation approximation 


Fig.7.2 Illustration of two approximations of a function g. In the left panel, we illustrate a piecewise 
constant approximation defined by the value of g at x = x; for the interval between x; to x;+1. In 
the right panel, we illustrate a piecewise linear approximation defined to coincide with g for x = x; 
and x = xj41. 


approximation is simply a constant function defined by the value of g at x = x;. The 
second is a linear function that coincides with g for x = x; and x = x;4. 

If we use these two approximations to estimate the integral of g from x = x; to 
X = Xj+1, We get 


А g(x)dx = Axg(x;) (7.14) 


and - i 
n g(x)dx « zA) + 2(Xi+1)), (7.15) 


respectively. Since we clearly have 


1 М-1 рх 
ах = ах, 7.16 
f вода 2J o(x)dx (7.16) 


we get two approximations of the integral from (7.15) and (7.16), respectively; 


1 M-1 
/ g(x)dx = Ax 2, g(xi), (7.17) 
апа | | "T 
n Уб + абын), (7.18) 


Here, the latter formula can be rewritten slightly to 


1 
f g(x) Ах Gn) + або) + g0) ++ кхм) 586). 119) 
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Table 7.2. Maximum errors of the Riemann sum approximation (Ев, (7.17)) and the Trapezoidal 


scheme approximation (Ет, (7.19)) to the integral h x?dx for some different values of Ax. 


Ax Eg Ер /Ах Ет Ет/Ах? 
0.1 0.0483 0.48 0.00167 0.17 
0.02 0.00993 0.5 6.67.10 0.17 
0.01 0.00498 0.5 1.67-1075 0.17 
0.002 0.000999 0.5 6.67-1077 0.17 
0.001 0.0005 0.5 1.67-10-7 0.17 


The approximation of the integral defined by (7.17) is referred to as a Riemann 
sum, whereas the approximation given by (7.19) is referred to as the Trapezoidal 
method of integration. In Table 7.2, we show the error when using these to schemes 


to approximate the integral 
1 
1 
Г х?ах = = 
0 3 


using several values of Ax. For the Riemann sum, the error is O(Ax), and for the 
Trapezoidal scheme, the error is O(Ax?). So, by using the same number of function 
evaluations, the accuracy of the integration scheme is improved from first to second 
order. 


7.2.3 Second Order Operator Splitting 


In order to introduce second order operator splitting, it is useful to introduce notations 
that we can use to simplify the formulations. We recall that the problem we want to 
solve is the following reaction-diffusion equation, 


ш = Uxx + f(u), (7.20) 
with initial and boundary conditions, 


u(0, x) = и0(х), (7.21) 
и(ї,0) = u(t, 1) = 0. (7.22) 


First, we let D(At) be an operator that evolves the solution of the diffusion equation 
one time step At ahead. If the function u = u(t,-) is known at time t, then u(t + At,-) = 
D(At)u(t, -) denotes the solution of 


Ит = Uxx, (7.23) 
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at time f + At, where u(t) is the initial condition at time t. Similarly, we let R play 
the same role for the reaction part of the equation. If u = u(t,-) is known, then 
u(t + At,-) = R(At)u(t,-) denotes the solution of 


и, = f(u), (7.24) 


at time f + At, where u(t) is the initial condition at time t. With this notation, we can 
rewrite the first order operator splitting derived above in a very compact manner. 
The step from f, to t,4; can be written 


и"! = (А) Р(Аг)и". (7.25) 
Ву simply repeating the process, we find that 
и" = (R(At)D(At))"u°. (7.26) 


Equipped with this notation, we can improve the accuracy of the operator splitting 
using the same approach as for numerical integration (7.19). In order to improve the 
accuracy of the operator splitting to second order, we now approximate one time 
step by 

u"*! = D(At/2)R(AtD(At/2)u". (7.27) 


By combining several steps, we get 
u” = [D(At/2)R(At)D(At/2) --- D(At/2)R(At)D(At/2)]v°. (7.28) 
Note that by the definition of D, it has the following useful property, 
D(At/2)D(At/2) = D(At), 


since applying D(At/2) twice simply means to solve the equation twice using the 
time step At/2 . By using this property, we can rewrite (7.28) as follows, 


u” = (D(A1/2)[R(At)D(AN) "RAND (Ar/2)) и. (7.29) 


So the only difference between the first order scheme (7.26) and the second order 
scheme (7.29) is that in the latter, we do a half time step in the first and last time 
steps. This is very similar to the half step in the first and last step of the numerical 
integration above; see (7.19) compared to (7.17). 


7.3 Operator Splitting Applied to a System of Reaction-Diffusion 
Equations 


In order to demonstrate the use of the operator splitting techniques introduced above, 
we revisit the FitzHugh-Nagumo system, 
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Uz = бох + соо — а)(1 — v) — cow, (7.30) 
ш, = b(v — dw). (7.31) 


We use the same constants as above, 
a = —0.12, сү = 0.175, cp = 0.03, b = 0.011, d = 0.55, 6 =5- 1077. (7.32) 


Here, the solution operator of the diffusion step 2) evolves 


0; = ÓUxx, (7.33) 
wr = 0, (7.34) 
and the R evolves 
ш = civ(v — a)(1 — v) — cow, (7.35) 
ш; = b(v — dw). (7.36) 


In Table 7.3 we show the error when the first order (see (7.26)) and second order (see 
(7.29)) algorithms. The solutions are computed by replacing D and Ж by the standard 
explicit schemes with a fine time step (At — 107^) and the errors are estimated by 
comparing the solutions to a numerical solution found using the fine time step and 
no operator splitting. As anticipated, the convergence of the two methods are first 
and second order. 


Table 7.3 Maximum errors of the first order (E, (7.26)) and the second order (E2, (7.29)) operator 
splitting techniques applied to the FitzHugh-Nagumo system (7.30)-(7.32). In each operator 
splitting step, the system is solved using standard explicit schemes with At = 107^ and Ax = 
0.01. The error is found by comparing the solutions to solutions found using a standard explicit 
scheme with At = 107^ and Ax = 0.01 without operator splitting. 


At Ei Е / At E» Ej[At? 
5 0.0205 0.0041 0.00612 0.00024 
2 0.00768 0.0038 0.0011 0.00028 
1 0.00384 0.0038 0.000296 0.00030 
0.5 0.00192 0.0038 748-107? 0.00030 
0.2 0.000765 0.0038 1.05-107? 0.00026 


7.4 Comments and Further Reading 


1. The example in Section 7.3 indicates the strength of operator splitting. If we 
have a good solver for the diffusion step, and a good solver for a system of 
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ordinary differential equations, these solvers can be combined to get a solution 
of the reaction diffusion problem. In this case, this is not really necessary, but 
there are extremely complex problems out there that are virtually impossible to 
solve without using operator splitting. 

2. Why does operator splitting work? This question is studied in some detailed in, 
e.g., [4, 5, 7, 8, 9, 12]. A detailed discussion of this topic is far outside of our 
scope. However, we can give a hint to why this works. To this end, we consider 
a linear system of ordinary differential equations 


и, = Au + Bu, 


where u is a vector and А and B are matrices compatible with the dimensions of 
u;they may represent discrete version of spatial derivatives as above. An explicit 
scheme for this system can be written on the form 


ит! = (I + AtA + АтВ)и", 
and using first order operator splitting, we get the scheme 


U"*1? = (I + ArA)U", (7.37) 
UM = (I + Amm. (7.38) 


By combining the two steps of the operator splitting scheme, we get 
U"*! = (I + AtA)(I + ALB)U". (7.39) 
This scheme can be expanded to read 
U"*! = (I + AtA + AtB)U” + At? ABU”. (7.40) 


Hence, in the step from f, to tn+1 the results of the two schemes differ by O(At?), 
and thus the difference summed over N ~ At^! steps is O(At). Therefore, the 
standard explicit scheme and the first order operator splitting scheme converges 
to the same solution as Af goes to zero. 

3. Thefirst order operator splitting is referred to as Godunov splitting [3] and second 
order splitting is referred to as Strang splitting [10]. In [2], these methods are 
used to split spatial parts of the differential equation. In that setting, the method 
is referred to as the method of fractional steps. 

4. In computational electrophysiology, operator splitting was introduced for the 
monodomain model in [6] and extended to the bidomain model in [11]. The 
error of the method was analyzed in [7], and higher order splitting methods were 
introduced in [1]. 

5. In this chapter, we have learned that operator splitting is an effective method for 
breaking down complex problems into simpler ones. This technique is widely 
used in mathematics and is particularly useful in computational mathematics, 
as it enables the reuse of code. However, the question remains whether a 


62 7 Operator Splitting 


coupled system is inherently more difficult to solve. It is possible to solve a 
reaction-diffusion equation using a fully coupled implicit scheme and Newton’s 
method, but when the problem involves large systems of equations defined on 
different domains and scales, and it can be impractical to solve the problem in 
a fully coupled manner. In such cases, operator splitting is often necessary to 
manage the complexity of the problem. 
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Part II 
Models of Electrophysiology 


A 
Check for 


Chapter 8 peas 
Membrane Models 


In the previous chapters, we introduced techniques that can be used to find numerical 
solutions of differential equations. The examples we considered were simple, 
theoretical differential equations without units. From this point forwards, we will 
look at how to apply the methods introduced in the previous chapters to differential 
equations that are set up to model aspects of electrophysiology. 

In this chapter, we will consider a type of model that is commonly used to model 
the dynamics across the membrane of excitable cells. We will start by introducing 
a model for the action potentials generated in neurons. More specifically, we will 
consider the Hodgkin-Huxley model for the action potential of the squid giant axon 
[10]. Then, we will introduce a similar model for a cardiac action potential — 
specifically, a model for the action potential of a rabbit ventricular cardiomyocyte 


[9]. 


8.1 The Hodgkin-Huxley Model 


The Hodgkin-Huxley model [10] consists of a system of four ordinary differential 
equations with the four unknowns v, m, h, and r !: 


x = —(Una + Ik +), (8.1) 
2" = 0.(1 — m) – Bum, (8.2) 
P an(1 ~ h) - Brh, (8.3) 
T =a,(1-r)- yr. (8.4) 


! The unknown function r is actually called n in the original formulation of the Hodgkin-Huxley 
model, but we will use the name г to avoid confusion with the index п used to denote the time step 
in the numerical schemes. 
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Here, the unknown function v has unit millivolts (mV) and represents the membrane 
potential, and Cm = 1 uF/cm? is a parameter representing the specific membrane 
capacitance. Furthermore, Ina, Jk, and / represent the current density through three 
types of ion channel: Na* channels, K* channels, and non-specific leak channels. 
The current densities are given in units of wA/cm? and are defined by 


Iva = guam? h(v — ома), (8.5) 
Ik = gkr (v — ок), (8.6) 
Ij = gy(v - uL), (8.7) 


where gna = 120 mS/cm?, gx = 36 mS/cm? and ет, = 0.3 mS/cm? are parameters 
representing the maximal conductance density of the different channel types, and 
UNa = 50 mV, ок = —77 mV and uj, = —54.4 mV are the equilibrium potentials of 
the channels. In addition, т?л and r^ represent the open probability of the Na* and 
the K* channels, respectively. The open probability of the leak channels is assumed 
to be 1 at all times. 

The unknown functions m, h and r take values between 0 and 1 and are governed 
by the equations (8.2)-(8.4). In these equations, @m, Bm, Œh, Bh, @-, and B, represent 


rates of the opening and closing of channel gates, are given in units of ms! and 
depend on the value of v. More specifically, they are defined by 
B yiv + уз) P -(v*ys)] yi 
Om = Т e-G*yDha* Bm = Уде ло, (8.8) 
= уле @%%)/» = Y10 
SACS , Bn = 1 + e—@md/n2’ (8.9) 
yixo + yia) - (o 
* p eomas Br = ушне "mom, (8.10) 
where the parameters yı—y18 are constants specified by 
yi 20.1 ms !mV^!, y2 = 40 mV, уз = 10 mV, (8.11) 
ya = 4ms"!, ys = 65 mV, ув = 18 mV, (8.12) 
ут = 0.07 ms !, ys = 65 mV, yo = 20 mV, (8.13) 
ую = 1 ms^!, yii = 35 mV, ур = 10 mV, (8.14) 
yis = 0.01 ms! mV, ум = 55 mV, yis = 10 mV, (8.15) 
у16 = 0.125 ms, ут = 65 mV, yis = 80 mV. (8.16) 


In our computations reported below, we will use the initial conditions 


v(0)=-60mV, т(0) = 0.1, А0) = 0.6, r(0203. (8.17) 
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8.1.1 An Explicit Numerical Scheme 


The differential equations of the Hodgkin-Huxley model can be solved numerically 
using the techniques we have applied for simple example equations in the previous 
chapters. More specifically, we can define a numerical scheme for the equations by 
replacing the derivatives in (8.1)—(8.4) by the standard difference, 


f(t + At) - FO 


9 8.18 
Р) А (8.18) 
and define an explicit scheme by 
Un+1 — Un 
Cn = —(INa(0n, Mn, hj) + Ig(on, rn) + LL(n)), (8.19) 
с 72 = Ain(Un)(1 — Mn) — Bim, (8.20) 
hy a hy 
== = ол(о„)(1 — hn) — Bi(vs)hy, (821) 
T = ay (on)(1 — rn) – В.(о), (8.22) 


where, as usual, vn, Mn, An, and rp are the numerical solutions at time t, = n X At. 
The scheme can be written in computational form as 


Unc] = Un — EUN NN t Ikn: rn) + Tn), (8.23) 
Mn+1 = Mn + At[ao (v, )(1 m mn) Е В.„(0һ)т], (8.24) 
hn+1 = hn + At[ag(os)(1 — hn) – Bn(Vn)hn], (8.25) 
Intl = ry + At[o (os (1 — r4) — Br (Un)rn]- (8.26) 


8.1.2 Error of the Numerical Solution 


In Table 8.1, we report the error, E,, of the numerical solution of v using the numerical 
scheme (8.23)-(8.26) for some different values of Ат. The error is computed by 
comparing the solutions at time f = 3 ms to the solutions found using a very fine 
time step (At = 1076 ms). As observed for simpler systems of equations in earlier 
chapters, we find that the error of the explicit scheme is close to proportional to the 
time step, Ат, applied in the numerical scheme. In other words, we have linear (or 
first order) convergence. 
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Table 8.1 Error of the numerical solution of the Hodgkin-Huxley model for different values of Ат. 
The error is defined as E, = |v — vy |, where v is the numerical solution at t = 3 ms for a very fine 
time step (At = 10-6 ms), and vy is the numerical solution at t = 3 ms for each of the values of At 
reported in the first column of the table. 


At (ms) E, (mV) Е, /At (mV/ms) 
0.01 0.982 98 
0.005 0.49 98 
0.001 0.0979 98 
0.0005 0.0489 98 
0.0001 0.0097 97 


8.1.3 Details of the Model Solution 


In Fig. 8.1, we show plots of the numerical solution of the Hodgkin-Huxley 
model computed using At = 0.001 ms. In the upper left panel, we have plotted 
v, representing the membrane potential. We see that an action potential is generated, 
starting with an increase in the value of v (depolarization) followed by a decrease in 
the value of v (repolarization), like in the simple FitzHugh-Nagumo model studied 
in Chapter 2. The action potential lasts for a couple of milliseconds. 

In the next three panels of Fig. 8.1, we show the value of the unitless m, ^, and r 
variables, and in the final three panels we show each of the current densities Дуа, Ik, 
and J. We see that Лу obtains negative values, while /k is positive. The current Д, 
obtains both positive and negative values, and these values are considerably smaller 
(in absolute value) than the values of Дуа and /к during the action potential. The sign 
of the current densities can be explained by revisiting the definitions, 


Іа = guam? (о = UNa), 
Ik = exri(v — UK). 
TL = 81(0 - yy), 
(see (8.5)-(8.7)) and recalling that gna, gx, gL, m, h and r are all positive. Thus, 


we see directly from the definition of the current densities that Jya is positive for 
v > vya = 50 mV and negative for v < 50 mV. Similarly, 7x is positive for v > vy = 


—T1 mV and negative for v < —77 mV, and Д. is positive for v > v = —54.4 mV 
and negative for v « —54.4 mV. 
Moreover, since Ср; = —(Iwa + Ik + 1), see (8.1), and Cm is positive, we deduce 


that a negative value of the sum (Дуа + Ik + JL) is needed for v, to be positive (i.e., 
for v to increase) and a positive value of the sum (Дуа + Jk + Д.) is needed for v; to be 
negative (i.e., for v to decrease). We can therefore conclude that in the beginning of 
the simulation, for v < —54.4 mV, both Jya and д, contribute to the depolarization 
of v, and that after v increases above —54.4 mV, Iya is solely responsible for the 
depolarization. However, as the value of v increases, the value of (v — vk) and r 
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Fig. 8.1 Numerical solution of the Hodgkin-Huxley model with initial conditions v(0) 2 —60 mV, 


m(0) = 0.1, A(0) = 0.6 and r(0) = 0.3. The system of equations is solved using the explicit scheme 
described in (8.23)-(8.26) with Ат = 0.001 ms. 


increases, which leads to an increased value of Ig. When /к + J, is more positive 
than Дуа is negative, v; becomes negative, and v starts to repolarize. 
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8.1.4 Upstroke Velocity and Action Potential Duration 


From the discussion in the previous subsection, it seems reasonable to expect 
that increasing а, e.g., by increasing the value of the parameter ема, will make 
the depolarization phase (upstroke) of the action potential more rapid and the 
repolarization slower (longer action potential duration)?. Similarly, we would expect 
that increasing /к, e.g., by increasing the value of the parameter gx, would make the 
repolarization of the action potential more rapid, and thus the duration of the action 
potential shorter. In Fig. 8.2, we show the upstroke and the action potential for a few 
choices of the parameters gna, gx, and гт. As expected, we observe that increasing 
£Na increases the upstroke velocity and the action potential duration, while increasing 
£k decreases the action potential duration. The adjustments of gg do not seem to 
have a significant effect on the computed upstroke or action potential duration. 


50 Upstroke 50 Action potential 
S 0 0 
E 
> 
-50 -50 
0.5 1 1.5 0 2 4 
50 50 
———D.75xg 
S o 0 K 
[3 1x Ik 
> 
-50 -50 ———1.25х9К 
0.5 1 1.5 0 2 4 
50 50 
——0.5х0 
© o 0 L 
E 1x9, 
> А 
-50 -50 ——15x9, 
0.5 1 1.5 0 2 4 
t (ms) t (ms) 


Fig. 8.2 Numerical solution of the Hodgkin-Huxley model for some adjustments of the parameters 
ма, BK, and gL. The applied adjustments are given in the legends оп the right-hand side of each 
row, and the remaining parameters are kept at their default values. To make the comparison easier, 
the timing is adjusted such that the time of the maximal upstroke velocity occurs at the same time 
for all the parameter choices. The system of equations is solved using the explicit scheme described 
in (8.23)-(8.26) with At = 0.001 ms. 


? See Section 2.3 (page 17) for definitions of the upstroke velocity and action potential duration. 
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8.2 A Parsimonious Model for the Action Potential of Rabbit 
Ventricular Cardiomyocytes 


In addition to the Hodgkin-Huxley model for neuronal action potentials, we will also 
consider a similar model for a cardiac action potential. More specifically, we consider 
a simple, so-called parsimonious, model for the action potential of rabbit ventricular 
cardiomyocytes from [9]. Note, however, that numerous alternative cardiac action 
potential models exist (see Section 8.3). The parsimonious rabbit model consists of 
a system of three ordinary differential equations with three unknowns, 


dv 

Cn = —Una + + Istim), (8.27) 
dm Mo-m 
a аа 2 
dt Tu л) 
dh ho-h 
Em . 2 
dt Th 99) 


Again, the unknown function v in units of millivolts (mV) represents the membrane 
potential, and Cm = 1 uF/cm? is a parameter representing the specific membrane 
capacitance. Furthermore, Дуа and /к (in uA/cm?) represent the current densities 
through Na* and K* channels and are given by 


Iva = guam h(v — ома), (8.30) 
Ix = gke P") (y — vy), (8.31) 


where gna = 11 mS/cm? and gx = 0.3 mS/cm? represent the maximal conductance 
densities of Na* and К? channels, respectively, and vya = 65 mV and vg = —83 mV 
are the equilibrium potentials of the two channels types. In addition, т?л and 
e P "7X represent the open probability of the Na* and K* channels, respectively. 
The parameter b has the value b = 0.047 mV~!. As in the Hodgkin-Huxley model, 
the unknown functions m and л take values between 0 and 1 and are governed by 
(8.28)-(8.29)? where 


1 
Te ттт Tm = 0.12 ms, (8.32) 
_ 1 _ 210 gon (v En) 
fis e EUER Th = [ЕЕ УК? (8.33) 


апа 


3 Note that the formulation of the equation for m used here (m; = (mœ — m)/Tm) corresponds 
to the formulation used in the Hodgkin-Huxley model (m; = a,,(1 — m) – mm) if we define 


1 _ and Mmo = ae The equation for h can also be rewritten analogously. 
m+Bm 
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Fig.8.3 Values of mæ, ho and e P G-v9) in the parsimonious ventricular rabbit model (8.27)-(8.35) 
as functions of v. 


Em = —41 mV, km = —4.0 mV, (8.34) 
En =-74.9mV, К, =4.4 mV, t? =6.8 ms, 6, = 0.8. (8.35) 


In addition, Istim is a stimulus current density given in units of uAlcm?. This stimulus 
current density is introduced because the initial conditions (see below) yield a system 
at rest, and the stimulus current is needed to trigger an action potential. The stimulus 
current is described in more detail in Section 8.2.1 below. 

In our computations, we will use the initial conditions 


v(0) = -83 mV, m(0) = 0, h(0) = 0.9. (8.36) 


8.2.1 The Stimulus Current Density 


In the model (8.27)-(8.29), Istim is a stimulus current density that is used to initiate 

an action potential by increasing the membrane potential enough to activate Jya. 

Specifically, the stimulus current density is given by 

asim, if t > tstim and f < иш + stim, 

Бит = . (8.37) 
0, otherwise. 


In other words, the stimulus current density is only nonzero in the period from tstim 
tO fstim + dstim. In our computations, we use astim = —25 uAlcm?, tstim = 50 ms, and 
dstim = 2 ms, unless otherwise specified. This turns out to be a sufficiently strong 
Istim to activate yg and thus initiate an action potential. 

As mentioned above, without an included stimulus current density, the membrane 
potential, v, is at rest (i.e., does not change with time) at the initial conditions, v(0) = 
—83 mV, m(0) = 0, А(0) = 0.9. This is because both /к and Jya on the right-hand side 
of (8.27) are equal to or close to zero. The current density Ig = gke PP) (y — uk) 
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is zero because v = —83 mV, which is equal to the equilibrium potential of the K* 
channels, vg, so (v — vk) = 0. Furthermore, the current density Дуа = guam? h(v — UNa) 
is zero because m is zero. From Fig. 8.3, we also see that m..(v) ~ 0 for v = —83 mV. 
Therefore, mo; = m, which means that d - e ж Оапа mis expected to maintain 
a value close to zero. Since the membrane potential is at rest at v = —83 mV, it can 
be referred to as the resting potential of the model. 

However, if we apply a negative Itim for some time, this will allow v to increase a 
bit. If we for example increase v to about —40 mV, then m..(v) increases to about 0.5 
(see Fig. 8.3), and чи = = becomes positive (if m ~ 0), leading to an increased 
value of m. For m > 0, we get a nonzero, negative Дуа, that will last as long as h > 0 
and (v — ома) > 0 (see (8.30)), and this negative Jya creates the upstroke of the action 
potential. 


8.2.2 An Explicit Numerical Scheme 


An explicit numerical scheme for the parsimonious ventricular rabbit model can be 
defined in almost exactly the same manner as in Section 8.1.1 for the Hodgkin-Huxley 
model by replacing the derivatives іп (8.27)-(8.29) by the standard difference (8.18). 
The computational form for the scheme reads: 


At 
Un+1 = Un — C. Uam. Mns hn) + ШАШУ, T Tstim(tn)], (8.38) 
i en a Оа (8.39) 
Tm(Un) 
Noo п) - hn 
hn+1 = hn + Ag re) = n (8.40) 
Тһ(Оп) 


8.2.3 Numerical Computations 


In Table 8.2, we report the error, E,, ofthe numerical solution of v using the numerical 
scheme (8.38)-(8.40) for some different values of Ат. The error is computed by 
comparing the solutions at time t = 10 ms to the solutions found using a very fine 
time step (At = 107? ms). As observed for the Hodgkin-Huxley model, we find that 
the error of the explicit scheme is close to proportional to the time step, Ar, applied 
in the numerical scheme. So again we have linear (or first order) convergence. 

In Fig. 8.4, we show the numerical solution of the parsimonious ventricular rabbit 
model solved using At — 0.001 ms. In the upper left plot we show the membrane 
potential, v, and we see that like in Fig. 8.1 for the Hodgkin-Huxley model, an 
action potential is generated. The action potential starts around the time when Istim 
is applied (at fstim = 50 ms) and seems to last for about 250 ms. This is significantly 
longer than the about 2-3 ms long neuronal action potential in the Hodgkin-Huxley 
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Table 8.2 Error of the numerical solution of the parsimonious ventricular rabbit model for different 
values of At. The error is defined as E, = |v — vy |, where v is the numerical solution at £ = 10 ms 
for a very fine time step (At = 107? ms), and vy is the numerical solution at t = 10 ms for each 
of the values of Ат reported in the first column of the table. In these computations we have used 
faim = 0 ms and dgim = 2 ms. 


At (ms) E, (mV) Е, [А (mV/ms) 
0.01 0.662 66 
0.005 0.322 64 
0.002 0.127 63 
0.001 0.0627 63 
0.0005 0.0309 62 
v m h 
50 1 1 
as 0 1 
> 
E E 0.5 = 0.5 
> .50 J 
-100 0 0 
0 200 400 0 200 400 0 200 400 
Iya Na Ik 
0 0 3 
t E = 
5 -100 1 5 -100 Б 2 
E 3 3 
w -200 | 3-200 RA 
_2 par — 
-300 -300 0 
0 200 400 51 52 53 0 200 400 
t (ms) t (ms) t (ms) 


Fig. 8.4 Numerical solution of the parsimonious ventricular rabbit model. The system of equations 
is solved using the explicit scheme described in (8.38)-(8.40) with At = 0.001 ms. Note that Jya is 
shown in two panels, one with the same time scale as the remaining panels, and one with the time 
scale zoomed in on the time of the peak current. 


model for a neuronal action potential. The next two upper panels show how the value 
of m and Л changes with time, and the plots in the lower panels show the two ion 
channel current densities, Jy, and Jk. The lower left panel shows Дуа in the same 
time scale as in the remaining panels, and the lower middle panel shows Jya zoomed 
in on the points in time when the peak current occurs. The lower right panel shows 
Ik. 
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8.2.4 Upstroke Velocity and Action Potential Duration 


In Fig. 8.5, we investigate how the upstroke velocity and action potential duration 
are affected by adjusting the parameters ома and gx in the parsimonious ventricular 
rabbit model. In the upper panel, we observe that decreasing gya leads to a slower 
upstroke and a shorter action potential duration, whereas a lower value of gx leads 
to a longer action potential duration. These observations are consistent with what we 
observed for the Hodgkin-Huxley model in Fig. 8.5. We also observe that the onset 
of the rapid upstroke happens slightly faster after the stimulus current is applied 
when gx is decreased. 


Upstroke Action potential 
50 50 
—— 0.5x 9ма 
= 0 — 0.75х 94a 
Е 0 1х ма 
> -50 — 1.25xg,, 
—— 1.5x 9ма 
-50 -100 
51.5 52 52.5 53 0 200 400 
50 50 
—= 0.5x 9% 
S 0 — 0.75xg, 
Е 0 1х9К 
> -50 —— 1.25x 9к 
——— 1.5x 9k 
-50 -100 
51.5 52 52.5 53 0 200 400 
t (ms) t (ms) 


Fig.8.5 Numerical solution of v in the parsimonious ventricular rabbit model for some adjustments 
of the parameters gna and gx. The applied adjustments are given in the legends on the right-hand 
side of each row, and the remaining parameters are kept at their default values.. The system of 
equations are solved using the explicit scheme described in (8.38)-(8.40) with Ат = 0.001 ms. 


8.3 Comments and Further Reading 


1. Following the publication of the Hodgkin-Huxley model [10] in 1952, many 
similar models representing the action potential of different cell types were 
published, building on the same basic principles. One example is the ventricular 
rabbit model described in Section 8.2. More recent models are often more 
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complex than the two membrane models considered in this chapter. For instance, 
a large number of additional membrane currents and dynamic intracellular ion 
concentrations are often included in the equations. Nevertheless, the resulting 
system of differential equations can be solved in exactly the same manner as 
seen for the Hodgkin-Huxley model in Section 8.1.1 and for the parsimonious 
ventricular rabbit model in Section 8.2.2. 


. А nice overview of the evolution of mathematical membrane models for 


cardiomyocytes is found in [1]. These models include models for different 
cell types, like Purkinje fibres (e.g., [4, 14, 19]), ventricular cardiomyocytes 
(e.g., [7, 16, 20]), and atrial cardiomyocytes (e.g., [3, 8, 15]). They are set up to 
represent cells from different species, like mice (e.g., [2]), guinea pigs (e.g., [6]), 
rabbits (e.g., [18]), and humans (e.g., [16]). Membrane models representing stem 
cell derived cardiomyocytes have also been introduced (e.g., [11, 12, 13, 17]). 


. The membrane models in the form we have considered in this chapter are 


expressed as a system of ordinary differential equations. The membrane 
potential, v, obtains a single value at each time step and is governed by a number 
of current densities. There is no spatial variation present in the model. In reality, 
the currents across the membrane happen through ion channels at different 
locations of the membrane. The current densities in the models in the form 
considered in this chapter could therefore be interpreted as the average current 
densities over a given area of membrane, for example, over the membrane of one 
cell. In this case, v could be interpreted as the averaged membrane potential of 
the cell. In the next chapters, we will combine the membrane models considered 
in this chapter with spatial models of electrophysiology. 


. The formulation of the Hodgkin-Huxley model given in Section 8.1 is adjusted 


for the membrane potential to have a resting state at v = —65 mV, and is taken 
from [5]. 


. Compared to the original publication of the parsimonious ventricular rabbit 


model [9], we have for simplicity used values of Т, and бр rounded off to one 
decimal point. 


. A cardiomyocyte action potential typically lasts several hundred milliseconds. 


In the parsimonious model considered above, it typically lasts between 200 ms 
and 500 ms (see Fig. 8.5). In Fig. 8.4, we note that the gradients of the solution 
are very sharp, and this indicates that we need to apply small time steps to pick 
up the main features of the solutions. From Table 8.2, we note that the error is 
about E, = 63 x At mV/ms. Thus, if we accept an error of about 1 mV, we need 
the time step to be less that 1/63 ms. If the action potential lasts for 500 ms, we 
need to perform at least 31,500 time steps to achieve sufficient accuracy. 


. The action potential of a neuron lasts only a few milliseconds, while that of a 


cardiomyocyte lasts several hundred milliseconds. The difference in duration is 
due to a difference in the density of Na* and K* channels in the membrane of the 
two cell types. In the models considered in this chapter, the neuronal model has 
a 10 times higher density of Na* channels and around a 100 times higher density 
of K* channels than the cardiac model. The high density of channels leads to a 
brief and clear signal transmission in the neuronal model. On the other hand, the 


References TI 


longer action potential of the cardiomyocyte is crucial in inducing an increase in 
intracellular Ca?* concentration (not part of the model above), resulting in the 
mechanical contraction of the myocyte and effective pumping of blood. Thus, 
the longer action potential is essential for optimal cardiac function. 
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Chapter 9 
The Cable Equation 


In Chapter 6, we studied a simple version of the cable equation, where a diffusion 
term was added to the FitzHugh-Nagumo equations. In this chapter, we will revisit 
the cable equation and go through a simple derivation of the model. In addition, we 
will consider the numerical solution of the cable equation for a neuronal axon with 
membrane dynamics modeled by the Hodgkin-Huxley model. 


9.1 Derivation of the Cable Equation 


In order to get a sense of the origin of the terms in the cable equation, we will here 
consider a simple derivation of the model. Similar derivations are found in, e.g., 
[1, 2, 4, 7, 9]. The derivation is based on dividing a cell (e.g., an axon) into a number 
of compartments in the x-direction, as illustrated in Fig. 9.1. 


— Ах 


Fig. 9.1 Left: Illustration of a cell separated into a number of compartments of length Ax. Right: 
Illustration of one of the compartments, j. The cross-sectional area in the x-direction is denoted by 
Ax, and the total membrane area is denoted by Am. The current across the membrane is denoted 
by Im, the current from compartment j to compartment j + 1 is denoted by /, and the current from 
compartment j to compartment j — 1 is denoted by /... 
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The left panel of Fig. 9.1 illustrates a cell divided into a number of compartments 
of length Ax, and the right panel illustrates an arbitrary inner! compartment number 
j. We will derive the cable equation by considering each of the electric currents 
flowing out of this compartment. 


Currents Between Compartments 


We assume that the currents that flow between compartments are governed by Ohm's 
law in the sense that the current from compartment j to compartment j + 1 is given 
by 
Ui j — Мї,]+1 
R , 
where u; у is the electrical potential the center of compartment j and и; j+1 is the 
electrical potential in the center of compartment j + 1, both in units of millivolts 
(mV). The subscript i is used to specify that we are considering the intracellular 
potential of the cell. Furthermore, R is the intracellular resistance between the two 
compartment centers in units of kilo-Ohm’s (kQ). This resistance can be expressed 
as (see, e.g., [5]) 


I, = (9.1) 


Ах 
R = —, 9.2 
TY (9.2) 
where Ax is the distance between the centers of compartments j and j + 1 (in cm), 
o; is the intracellular conductivity (in mS/cm), and A, (in cm?) is the cross sectional 
area of the cell (see Fig. 9.1). Inserting (9.2) into (9.1), we get 


Ui, j 7 Mijn 


L = oj; Ax Ax , 


(9.3) 
and, following the same steps, we get that the current from compartment j to 
compartment j — 1 is given by 


Ui j — Ui j-1 


І = оАх 
Ах 


(9.4) 


The unit of I, and /_ is micro-Amperes (uA). 


Membrane Currents 


We assume that the membrane acts as a capacitor that can store charge and that 
this property can be modeled like in the membrane models in Chapter 8. More 
specifically, we assume that the total current across the membrane of compartment 


! Here, an arbitrary inner compartments refers to any of the compartments except for the rightmost 
or leftmost compartments. 
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j is given by? 


Ov; 
In = Am Спот + Поп |. (9.5) 


Неге, Am is ће membrane area of compartment j (in cm?), given by 
Am = NmAx, (9.6) 


where nn is the circumference of the compartment (in cm) and Ax is the length of 
the compartment (in cm). Furthermore, Cm is the specific membrane capacitance 
(given in units of uF/cm?), and v; is the membrane potential (in mV) defined as the 
potential difference 

U = Uj — Ue, (9.7) 


where и; is the intracellular potential in the compartment and ue is the extracellular 
potential outside of the compartment. The term Cm Au is called the capacitive current 
density and represents the current density to and from the collection of charges stored 
by the membrane capacitor. Moreover, [jon is the current density (in uA/cm?) through 
ion channels in the cell membrane. For example, if we assume that the membrane 
dynamics are modeled by the Hodgkin-Huxley model (see Section 8.1), this current 
density is given by 

Lon = Ina + Ig + L. (9.8) 


Sum of Currents 


In order to derive the cable equation, we assume that Kirchhoff's current law applies 
in each compartment. In other words, we assume that the sum of all of the currents 
out of the compartment is zero. This gives 


I +1 +1 = 0, (9.9) 
and inserting (9.3)-(9.6), this yields 


Ui j — Ui, jl 


ve pud u— 
LEE Ax Ax 


Бас di 2 =0, (9.10) 


which can be rearranged to 


до; oj Ax Ui j-1 = 2ui,j + Ui jel 
E о C J = с. (9.11) 


2 [n (9.5) we recognize the terms C,,,v; and Tion (see (9.8)) from the membrane models considered 
in Chapter 8. In those models, we ignored all spatial variation and assumed that /,4 was the only 
current into or out of the cell. In order for the sum of the currents to be zero, we therefore ended up 
with models on the form Ср; = —lion. 
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Now, we insert an assumption that the extracellular potential is zero everywhere such 
that v = u; (see (9.7))?. In that case, the model reads 


до; 0; — 20; + Viy 
Са = 8 һө (9.12) 
where А 
QR. (9.13) 
Nm 


From Chapter 3 (see page 23), we recall that 


v(t, x — Ax) — w(t, x) + v(t, x + Ax) 


pw; (9.14) 


Uxx(t,X) © 
for a sufficiently small Ax. We therefore assume that Ax in (9.12) is very small, and 
obtain the cable equation 


Cm = f- — lion: (9.15) 
x 


Since the considered compartment j was chosen as an arbitrary inner compartment, 
we conclude that the equation (9.15) applies everywhere along the cell, except at 
the left and right boundaries. The boundary conditions are considered in the next 
subsection. 


9.1.1 Boundary Conditions 


Intuitively, for the leftmost compartment of the cell (see Fig. 9.1), the current from 
this compartment to the non-existing compartment to the left, /_, should be assumed 
to be zero. In other words, 


Ui j — Ui j-1 


І = оАх 
" Ax 


= 0. (9.16) 


From previous chapters (see, e.g., (1.6) in Chapter 1), we know that 


ди uj(t, x) — uj(t, x — Ax) 
ES 17 
дх з) Ах ко 


for a sufficiently small Ax. Assuming that Ax is very small, (9.16) can therefore be 
translated to the boundary condition 


ди; 
о, (9.18) 
дх 


3 Note that the assumption that the extracellular potential is zero could be replaced by alternative 
assumptions (see the Section 9.4 for more details). 
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where we have assumed that the leftmost boundary of the cell is located at x = 0. 
Dividing by o;A;. on both sides of the equation and using the assumption that the 
extracellular potential is zero, which gives u; = v, we get the boundary condition 


22 (0) = 0. (9.19) 


Inserting this into the discrete version of the cable equation for the left boundary, we 
get 
Ov Í —U| +02 
ðt др 
A similar argument for the right boundary of the cell, at x = L, gives the boundary 
condition 


= Tion. (9.20) 


д 
S Lys (9.21) 
дх 
and the discrete equation 
Ovum —UM + UM-1 
Cy arum 6 AX = Tion- (9.22) 


9.1.2 Geometry 


The value of 6 (see (9.13)) in the cable equation depends on the considered geometry. 
If we consider a rectangular cuboid, as illustrated in Fig. 9.1, with width w in the y- 
and z-directions, we have A, — u? and Nm = 4ш, which gives 

OjÁ, | шо; 


=— = ——. 2 
i Пт 4 и 


Similarly, if we consider a cylinder with radius r, we have A, = ar? and nm = 2л, 
which gives 
ojAx то 


ô= ~ = : 9.24 
"m ? (9.24) 


9.1.3 Additional State Variables 


As mentioned in the derivation above, the term Jion represents the current density 
through ion channels in the cell membrane. This current density could, for example, 
be modeled by the Hodgkin-Huxley model: 


lion = Jna + Ik + Д, (9.25) 
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(see Section 8.1). Here, the formulation of Jy, and /к involve the additional state 
variables m, h and r, which are governed by (8.2)-(8.4). In order to define the current 
density, Доп, we therefore need to include these equations in the model, where m, h 
and r are functions of both ѓ and x. We thus get a system of equations of the form 


Ov д?р 
mar = 6—5 — lion, 2 
© дї 259 vee) 
P ант) - f (9.27) 
h 
оп = ay(1 — h) – fih, (9.28) 
x =a,(1-r)- В,г. (9.29) 


9.2 Numerical Schemes 


We will consider two alternative numerical schemes for the cable equation. First, 

a straightforward explicit scheme and then an operator splitting scheme, treating 

the diffusion part of the system implicitly. In both schemes, we seek numerical 

approximations to the solution in the M spatial points x; — (j — 1) x Ax, for 
L 


j = Lh... M, at the time points tn = n x At for n = 1,...,N. Here, Ax = ўт and 


At = T where L is the length of the domain and T is the total simulation time. 


9.2.1 An Explicit Numerical Scheme for the Cable Equation 


We consider a numerical scheme for the cable equation that is based on this discrete 
spatial version of the equation considered in the derivation of the model, i.e., (9.12). 
To define а numerical scheme, we have to replace the remaining derivative бр bya 
difference, and we choose the usual difference (see, e.g., (1.5) on page 4). In addition, 
similar replacements of derivatives by differences are inserted into (9.27)-(9.29), and 


we get the explicit scheme 


vit Е vy ve Е 2v7 $ Vi n „п pn „n 
C x = $ re — Tron; ‚т, hit; ), (9.30) 
m”+! — m” 
a = mU = mj) - Bw mj, (9.31) 
йы 
2: = ол(оу X1 = hj) — Bi oj hj, (9.32) 
p" - г? 
CMS = arv; A — rj) - В.(о; )r; (9.33) 
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for j = 2,...,M — 1. For j = 1 and j = M, we replace the right-hand side of (9.30) by 
the right-hand sides of (9.20) and (9.22), respectively. This scheme can be rewritten 
to computational form 


where 


At At 
ME 0 "Е | ib c, o m ry 
met! = т" + At[æm (v (1 — т") — Bo" )m" |, 
pn =h" + Atlan(v")\(1 _ h”) — By (v")h"], 


rl = т Ацан (07) — r”) - Be(o")r" 


the matrix А is defined by 


1 Mi hī 
v! m” h” 
=| | m=| |, ata] ? |, r= 

UM mw hy 

-11 0. 0 

1-21 0 0 

A-.9|01-21 0 0 

Ax? " | 
0 0 1-1 


and J is the M x M identity matrix. 


9.2.2 An Operator Splitting Scheme for the Cable Equation 


(9.34) 


(9.35) 
(9.36) 
(9.37) 


(9.38) 


(9.39) 


A potential disadvantage of the explicit scheme defined above is instability issues, 
like observed in Chapter 3 (see Section 3.4.3). In an attempt to avoid these issues, 
we therefore also define an operator splitting scheme for the cable equation, with an 
implicit treatment of the diffusion part of the system. In other words, we divide the 
system into two parts for every time step by first solving 


Cin 


Ov _ 50” дт ðh дг 
ðt — Ax?’ дї , or! Ot 


using an implicit scheme and then solving 


(9.40) 
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Cm = = —lion. (9.41) 
one аһ(1 = m) – Ват, (9.42) 
дї 

һ 
оп = ол(1—Һ)— Bh, (9.43) 
x =a,(1-r)- brr, (9.44) 


using an explicit scheme. In computational form, using the vectors Uy, Mn, hy, Ty as 
defined in (9.38) and the matrix A as defined in (9.39), the scheme for the first step 
reads 


[r- aeu" m" eg, АЗ шд", үт? шүл, (945) 


т 


and the scheme for the second step reads 


pH = y”+!/2 Б A liu, тен, pers vnd (9.46) 
m" E m”+!/2 + Atlam(v"*/)(1 _ m”+!/2) _ Bs (ot yn), (9.47) 
pn = p? d At [ay (o *  2)(1 = А) =p EP (9.48) 
pn - pn + А о, (о?у 2 r”+1/2) pag uen. (9.49) 


9.3 Numerical Computations 


We apply the two numerical schemes to a simple test case for the cable equation 
with membrane dynamics modeled by the Hodgkin-Huxley model. We consider a cell 
shaped as a rectangular cuboid with length L = 0.5 cm and width ш = 0.001 cm, and 
we assume that c; = 4 mS/cm. This gives 6 = 0.001 mS (see (9.13)). Furthermore, 
the initial conditions are given by 


—50 mV, f = 0.05 em, 

(0, x) = mV, for x < cm (9.50) 
—65 mV, for x » 0.05 cm, 

т(0, х) = 0.1, h(0,x) = 0.6, r(0,x)-2 0.3. (9.51) 


9.3.1 Stability 


In Fig. 9.2, we show the numerical solution of the problem found using the explicit 
scheme with Ax = 0.001 cm and Ar = 0.001 ms at some different points in time. We 
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Fig. 9.2 Numerical solution 
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model at five different points -100 | | | 1 
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0 0.02 0.04 0.06 0.08 0.1 
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0 0.02 0.04 0.06 0.08 0.1 
x (cm) 


observe that already at the first time steps, we get unreasonable oscillations close to 
the discontinuity in the initial conditions (at x = 0.05 cm, see (9.50)). 

In Fig. 9.3, we have solved the system of equations using the explicit scheme with 
a smaller value of At (At = 0.0002 ms). In this case, we avoid the oscillations, and we 
seem to get a reasonable traveling wave solution. However, because of the small time 
step, the numerical computations are quite slow. For this reason, we also try to solve 
the system using the operator splitting scheme described in Section 9.2.2. Using this 
scheme, we are able to use a time step of At = 0.02 ms, and still get a solution that 
is very similar to the one that we got using the explicit scheme with a fine time step 
(compare the solid blue and dotted orange lines in Fig. 9.3). Even though we have to 
solve a linear system of equations in order to find the solution in the first step of the 
operator splitting scheme, the algorithm is able to find the solution much faster than 
the explicit scheme because of the large difference in the required time step. More 
specifically, for the laptop computer that we used to perform the computations, the 
operator splitting scheme was about 25 times faster than the explicit scheme. 
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Fig. 9.3 Numerical solution Бб #= 0 тѕ 
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9.3.2 Conduction Velocity 


In Chapter 6, we considered a unitless version of the cable equation with membrane 
dynamics modeled by the FitzHugh-Nagumo model and observed how the conduction 
velocity, CV (i.e., the velocity with which the traveling wave moved though the 
domain) depended on two of the model parameters. In Fig. 9.4, we show the results 
of a similar experiment for the cable equation with membrane dynamics modeled by 
the Hodgkin-Huxley model. We vary the value of the three parameters gna, gx, and 
gL representing the maximal conductance density of the three membrane currents 
of the Hodgkin-Huxley model (see (8.5)-(8.7)). In addition, we vary the value of 
the cell width, w, and the intracellular conductivity, o, used to define 6 in the cable 
equation (see (9.23)). 
We define the conduction velocity as 
0 = ХІ 


У = ——— .52 
cv = 2л, (9.52) 
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Fig. 9.4 Conduction velocity (CV) computed from the numerical solution of the cable equation as 
described in (9.52) for a few adjustments of parameters. The parameter values not explicitly stated 
on the x-axis are set to their default values. The numerical solution is found using the operator 
splitting scheme described in Section 9.2.2, with Ax = 0.001 cm and At = 0.02 ms. 


where x; = 0.2 cm and x» = 0.4. Furthermore, tı and f» are the points in time when 
the value of v first increases to a value v > 0 mV, in the spatial points x; and x», 
respectively. 

In Fig. 9.4, we observe that increasing the value of gna, с; or the cell width, 
w, significantly increases the conduction velocity. Increasing gr, also increases the 
conduction velocity somewhat, whereas increasing gx decreases the conduction 
velocity. 


9.4 Comments and Further Reading 


1. In the derivation of the cable equation given in this chapter, we assume that 
the extracellular potential is zero. However, if we instead assume that the 
extracellular potential is another constant, C, different from zero, we end up 
with the exact same formulation of the cable equation. In that case, we would 
have и! = рі + С (see (9.7)), and inserting this into (9.11), we get 


дрі Е vi! e C -2( + С) cv +С 
Cm ðt 2 Ax? ES 
рі – 2p) + vi 
= ó—AÁd — — Доп. (9.53) 


90 
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which is exactly the same as (9.12). 


. The cable equation can also be derived by assuming that the value of the 


extracellular potential can vary in the x-direction, but not in the y- and 
z-directions. This leads to a different formulation of the constant 6 in the cable 
equation. See, e.g., [4] for more details on this version of the cable equation. 


. In Section 9.1.1, we derived boundary conditions for the cable equation by 


considering the leftmost and rightmost compartments of the cell. However, 
since it might be reasonable to assume that the left and right sides of the cell 
are covered by membrane (see Fig. 9.1), we might also wish to assume that the 
amount of membrane area is larger for these compartments. More specifically, 
for the leftmost and rightmost compartments, we get 


Ax 
Am = Ax + Am = Ax Fe + n») А (9.54) 
Ах 


Inserting this into the discrete version of the cable equation, we get the adjusted 
i Ax 
óp = rx (9.55) 
Ах + Пт 


for the left and right boundaries. In the code associated with these notes, we have 
included an example simulation where we compare the solution of the system 
using this adjusted 6, at the boundary to the case where the default 6 is used 
everywhere. By running that code, we see that the solutions for these two cases 
are indistinguishable. 


. In this chapter, we considered the cable equation for modeling the spread of 


an electrical signal along a neuronal axon. However, the cable equation (9.15) 
can also be used to study the spread of an electrical signal along a cardiac fibre 
composed of a row of cardiomyocytes connected to each other by gap junctions. 
In this case, the increased resistance for the currents across the gap junctions 
has to be taken into account, either in an averaged manner (see, e.g., [6]) or by 
having a ó that depends on x in a way such that the discretely located cell-to-cell 
connections are represented (see, e.g., [3]). 


. In(9.10) we derived a compartmental version of the cable model equations. Then, 


we passed to the limit in Ax and obtained the differential equation (9.14). Finally, 
we replaced the derivatives by differences, and obtained the scheme (9.30) which 
is, more or less, the compartmental model we started by deriving. So, was the 
differential equation just a detour? Do we need it? No, we don't really need it, 
but it is customary to phrase models in terms of partial differential equations. 
The main reason for this may be tradition, but there are practical reasons as 
well. When the problem is formulated as a partial differential equation, we have 
a large collection of methods that can be applied; finite difference methods, 
finite volume methods, finite element methods, boundary element methods, and 
many more. Also, the problem can be expressed in a very compact manner using 
differential equations whereas the compartmental form (or finite difference form 
for that matter) is clunky and much harder to read for complex systems involving 
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many equations. Finally, differential equations can be analyzed mathematically 
using tools that are harder to apply when the equations are discretized. But 
the major disadvantage of the formulation as a differential equation is that it 
is not straightforward to solve using computers — the form and popularity is 
inherited from a time where math was done using paper and pencils. An attempt 
to introduce partial differential equations using both continuous and discrete 
approaches can be found in [8]. 
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A 
Check for 


Chapter 10 гызы, 
Spatial Models of Cardiac Electrophysiology 


In Chapter 8, we introduced mathematical models of the action potential across a 
cell membrane. Understanding the properties of the cell membrane is absolutely 
essential in order to understand the electrophysiology of excitable cells. However, 
some essential properties can only be studied in spatially resolved models; i.e., in 
models representing spatial variation across a single cell or a collection of cells. 

Every heartbeat is based on an electrochemical wave traversing the whole cardiac 
muscle. Strong perturbations to these waves are referred to as arrhythmias. These 
disturbances can seriously disrupt the contraction of the heart muscle and are 
therefore very dangerous and potentially lethal. Cardiac fibrillation refers to a state 
of the heart where the contraction is completely unsynchronized, leading to severely 
reduced pumping functions. In mathematical modelling, fibrillation can only be 
studied in spatially resolved models, and hence the membrane models introduced 
above are inadequate. But these models are excellent building blocks in models 
representing collections of cells. 

A first step towards spatially resolved modeling was presented in Chapter 9, 
where we discussed the cable equation. This model is often used to study a strand 
of cardiac cells, but since the model is inherently one-dimensional, it has limited 
relevance for complex spatial phenomena like cardiac fibrillation. Here, we will 
present the celebrated bidomain model. It is considered to represent the gold standard 
of mathematical models of cardiac electrophysiology and dates back 50 years. From 
the bidomain model it is easy to derive the somewhat simpler monodomain model. 
This model is easier to deal with in terms of numerical solution and is therefore 
frequently used as a replacement of the more correct bidomain model. 

Here we will just present the models and show that the techniques we have 
derived above can be used to obtain numerical solutions of both the monodomain 
model and the bidomain model. In the notes below we will point to literature where 
the models are derived, and we will point to better numerical methods and to 
applications of the models. Again, we will just consider the finite difference method 
in order to keep things as simple as possible (but not simpler) and therefore we 
will use two-dimensional squares as computational domains. However, excellent 
open-source software tools are available for numerical simulations based on the 
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bidomain and monodomain models using the finite element method. Using the finite 
element method, more realistic geometries can be applied. 


10.1 First, the Diffusion Equation, Again 


For purely technical and notational reasons, we will first briefly consider the 
two-dimensional diffusion equation. This is helpful in order to understand the mesh, 
the matrices and the vectors that we use to solve the bidomain and monodomain 
equations below. 

Consider the following initial and boundary value problem (unitless), 


ди@,»у) _  Pult.xy) „дих, у) 


= 10.1 
дї дх? ду? SUD 


for Q : (x, y) € (0,1) x (0,1) and for < T. We let д0 denote the boundary of 
the computational domain, О. At the boundary, д0, we use the Neumann boundary 
conditions 

Ou(t,0,y) _ Ou(t, 1, y) Ou(t, x,0) Ou(t, x, 1) 


0, =0, ———=0, ——=0, 10.2 
Ox Ox ду ду (10.2) 


and, in addition, we define the initial condition 
u(0, x, y) = u° (x, у), (10.3) 


where w(x, y) is a given function. In (10.1), с denotes a strictly positive (given) 
constant. In order to derive a numerical scheme for this problem, we proceed as 
usual by replacing derivatives by difference. By using the difference approximations 
introduced in Chapter 3, we get the following scheme!, 


k, k, k-1, k, kl, k,j-1 k,j*l 
: Lag E EV (10.4) 
At Ax Ay 


utl oan. и" .—2и".+и",. ик _ү—2ир,+ир 


Here, иг denotes an approximation of u(tn, Xk, уу) where t, = nAt, xy = (k — Ах 
and y; = u — 1)Ay with Ах = 1/(M, — 1) and Ay = 1/(M; – 1) for sufficiently large 
integers M, and My. 

We noted above (ses page 43) that it is quite useful to formulate these numerical 
schemes using matrix/vector notation. However, here the unknowns are given by 
uy j and thus it is not straightforward to use such a notation. We therefore introduce 
a one-dimensional numbering of the two-dimensional problem. Specifically, we 
define the new index i = M,(j — 1) + k. Then, any two-dimensional vector z with 
components {zx,;} can be rewritten as a one-dimensional vector z; where i runs 
from 1 to M = M, x My. The components of the one-dimensional vector are given 


! [t may be useful to note the similarity with the one-dimensional case; see (3.13) at page 23. 
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Matrix column number 


by см, (j-1)+k Where j and k runs from 1 to М, and My, respectively. Using this 
numbering, we can rewrite the scheme (10.4) as follows, 


п+1 п n _ п п п = n n 
си O Upa 2U te Шм, 20 t+ Mie, 
= с Л +0 2 : (10.5) 
At Ax Ay 


и 


By defining px = с /Ax? and py =a/ Ay’, we can rewrite the scheme as follows, 


ae =u; + At [ох(и? 4 + их) — 2(px + pyu; + py(uj_y, + up.) . (10.6) 


Now, we can define an M x M matrix А where non-zero elements of a typical row i 
in the matrix are given by 


aj,i-M, = Py» (10.7) 
Qj,i-1 = рх, (10.8) 
аы = —2(рх + ру), (10.9) 
йл+1\ = px, (10.10) 
diieM, = Py. (10.11) 


When the index i corresponds to a boundary point or to a corner point in the mesh, 
there are exceptions. The detailed definition of all the elements of the matrix is found 
in the software associated these notes, and a figure showing the structure of the 
matrix is depicted in Fig. 10.1 in the case Ax = Ay = 0.2 and ø = 0.04, which gives 
Px = Py = 1, Mx = My = 6, and M = MxM; = 36. 

With this notation at hand, we can write the scheme (10.6) in the convenient form 
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ит! = (I + AtAyu”. (10.12) 


Furthermore, we can define the implicit scheme, 


= Аи"+!, (10.13) 
which can be rewritten as, 
(1—- AtAyu" = и", (10.14) 


where, as usual, / denotes the identity matrix. We can also define the second order 
midpoint scheme as follows, 


1 

XR UR А" + и"), (10.15) 
ог n A 
t t 

\! = $4 u^ = \ + А) и". (10.16) 


10.2 The Bidomain Model 


The bidomain model with membrane dynamics modeled by the parsimonious 
ventricular rabbit model [6] can be written in the form 


à 

Y Cn + Ton(v,m,h)) = V - (o; Vv) + V - (o; Vue), (10.17) 
0 2 V- (o; Vv) + У. ((о; + с„)УМи„), (10.18) 

dn ты-т 
Bu. oc 10.19 
dt Tu ( ) 

dh hs-h 
Sa 10.20 
dt Th ( ) 


We consider this system for Q : (x,y) є (0,L) х (0, L) and for t < T. The spatial 
coordinate is given centimeters (cm) and time is given in milliseconds (ms). The 
unknown functions v and ue are the membrane potential and the extracellular 
potential, respectively (in units of mV). Note that v = и; — ue, where u; is the 
intracellular potential. We can use any pair of variables among v, и; and и, as prime 
variables, but it is most common to use v and ue. The bidomain model can be derived 
by assuming that the membrane potential, v, the intracellular potential, и;, and the 
extracellular potential, ue, are all defined everywhere in the domain and following 
the same steps as those used to derive the cable equation (see Chapter 9) without the 
assumption of ue = 0 (see Comment 2 in Section 9.4). 

In the bidomain model, c; and с, denote the conductivities of the extracellular 
and intracellular spaces, respectively. The conductivities are tensors allowing the 
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Table 10.1 Parameter values used in the bidomain model simulations. 


Parameter Value 

Cu 1 nu F/cn? 
Ti x Oi,y 3 mS/cm 
Oe,x, Се,у 10 mS/cm 
X 2000 cm! 
Ly, Ly 1 ст 

Ах, Лу 0.025 ст 
At 0.01 ms 


conductivity to vary according to the spatial directions. Specifically, we have 


бз = Се,х 0 inus Сх 0 
AEN су)” PU Siy)’ 


In general, the conductivities can vary in space, but here we will assume that they аге 
constant. Furthermore, C,, is the specific membrane capacitance, and y denotes the 
surface-to-volume ratio of the cell membrane. The sum of the ion current densities 
across the membrane are given by 


lios (v, m, h) = Ina(v,m, h) + Ik(v) + Гей. (10.21) 


The specific formulations for the current densities JN, and /к are as specified in 
Chapter 8 (see page 71), and so are the associated gating functions те, hoo, Tm and 
тһ. The current density Istim represents a stimulus used to start the electrical wave. It 
is given by 


Айт, if 2 нт and t < fstim + dstim, 


Ttim(t, x, y) = and yx? + y? < lim, (10.22) 


0, otherwise, 


where а;іт = —25 uA/cm?, tsim = О ms, @ т = 2 ms, and Бл = 0.25 cm. The 
initial conditions are as specified for the parsimonious rabbit model in Chapter 8 
(see page 72), and, for simplicity, we use the boundary conditions 


Ue(t,0, y) = ue(t, Ly, y) = ue(t, x,0) = ue(t, x, Ly) = 0, (10.23) 
Ou;(t,0, y) -0, Oui(t, Ly, у) -0, ди; (1, х,0) 21 ди(ї, x, Ly) e 
Óx Ox Oy ду 


(10.24) 
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10.2.1 Operator Splitting for the Bidomain Model 


In order to solve the somewhat intimidating system (10.17)-(10.20) we use, more 
or less, all the tricks we have introduced above. The main technique, however, is to 
break the complex problem into parts that we are able to deal with. The first step 
along that path is to apply operator splitting. We start by assuming that the complete 
solution vector given by (v, uc, m, h) is known at time t = f, = n X At, and we want 
to compute the solution at time t = ї„+1. We do this in two steps. First we solve the 
following system of ordinary differential equations, 


д 

Ca 22 = —lion(v,m, h) (1025) 
дї 
ат То = Т 
ar 10.2 
dt Ta (10.26) 
dh hoy —h 
ша 10.27 
dt Th (10.27) 


By solving these equations with ¢ ranging from 1, to tn+1, we obtain a solution that 
we denote (v, ue, m, h)"*!?. In the next step, we solve the spatial part of the equation, 


Cn = V (a Vv) + V - (o; Vue), (10.28) 
02 V-(o;Vv) + V - (с; + Ge) Vue), (10.29) 
dm 
E 10. 
2i 0, (10.30) 
dh 
HE 10.31 
ET 0, (10.31) 


where the initial conditions at time t = f, is given by (v, ue, т, hy*V ? Here, we 
immediately note that m"*! = т"*!?? апа h”*! = n”*!/2. In order to compute v and 
Ue at tn+1, ме need to solve the linear system, 


2 = V-(a;Vv) + V - (o; Vue), (10.32) 


02 V-(o;Vv) + V · (о; + Ge) Vue), (10.33) 


where the initial condition for v is given by v”*!/?, 


10.2.2 Finite Difference Approximation 


The task at hand is now to solve the system of ordinary differential equations 
given by (10.25)-(10.27) and then solve the partial differential equations given by 
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(10.32)-(10.33). In order to do this we introduce a mesh as above. Specifically?, we let 
uy j and Ur j denote approximations of u(tn, Xk, уу) and v(t,, Xk, уу) where t, = nAt, 
хк = (k — 1)Ax and y; = (j – 1)Ay with Ax = L,/(M, — 1), Ay = L,/(M, - 1) 
and At = T/M, for sufficiently large integers My, Му and M;. Again, we use the 
mapping from a two-dimensional notation to a one-dimensional vector notation by 


defining i = M,(j — 1) + k. 


Explicit ODE Step 


By using this notation, we can write an explicit scheme? for solving the ordinary 
differential equations given by (10.25)-(10.27) as follows, 


1 1 


At 
or = of - Hs fs mp hz) (00) + аа (охо у) 004) 


Moo(vj') — m? 


m? 2 = m; + Oe ` (10.35) 
ho(v?*) — n" 
peel? = т Ar M (10.36) 


where 1 <i < M, X My, and n = 0,..., M; – 1. 


Implicit PDE Step 


In the implicit PDE step we first note that m"*! = m"*!/? and h"*! = h”+!/2, For the 
systems of PDEs given by (10.32)-(10.33) we use the same discretization as we used 
in (10.4) for the diffusion equation. This leads to the definition of two matrices, A; 
and Ae, with typical elements given by (10.7)-(10.11). For Ae, the typical elements 
of the matrix are defined by using px = с, / Ax? and py = с„/Лу? in (10.7)-(10.11). 
Similarly, the typical elements of A; are defined by using px = o;/Ax? and ру = 
cil Ay? in (10.7)—(10.11). With this notation, we are ready to define the following 
scheme for the spatial part of the bidomain equations, 


рп! _ р"+1/? 
At 


0 = Аџи"! + (А; + Aul. (10.38) 


= А "+! + Ayu"! (10.37) 


This can be rewritten as a block matrix-vector system as follows, 


? We replace ue by u in order to reduce the load of subscripts. 


3 This is exactly the same scheme as we used for the parsimonious rabbit model in Section 8.2, only 
that we here need to solve the ODE system in all the computational nodes. 
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[= At Aj EN Ai pU prti/2 
Cn i XC КД À [ (10.39) 


10.2.3 Traveling Wave Solution of the Bidomain Model 


In Fig. 10.2, we show the solution of the numerical scheme for the bidomain model 
described above, with parameter values specified in Table 10.1 and membrane 
dynamics modeled by the parsimonious rabbit model described in Chapter 8. The 
upper panel shows the membrane potential, v, at four different points in time, and 
the lower panel shows the extracellular potential, ue, at the same time points. We 
observe a traveling wave solution initiated by the stimulus current applied in the 
lower left corner of the domain, leading to an increased value of v. This increase 
in v gradually spreads through the domain, and at t = 20 ms the wave has traveled 
almost all the way to the opposite corner of the domain. 

From this traveling wave, we can compute the conduction velocity, for example, 
by 


D — 2 
CV = ү х) ton yi) l (10.40) 
2-t 


where x, = у = 0.4 cm and x2 = ух = 0.8 cm. Furthermore, tı and t are the points 
in time when the value of v first increases to a value v > —20 mV, in the spatial points 
(ху, y1) and (x2, y2), respectively. In that case, we find that the conduction velocity 
in this default case is about 54 cm/s. In Fig. 10.3, we investigate how the conduction 
velocity depends on the conductivity parameters, c; and cz, and we find that the 
conduction velocity increases if oj or ст, are increased. 


х (cm) 


Fig. 10.2 Solution of the bidomain model at four different points in time. The upper panel shows 
the membrane potential, v, and the lower panel shows the extracellular potential, ue. A traveling 
wave solution is initiated by stimulating the cells in the lower left corner. 
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Fig. 10.3 Conduction velocity computed from the solution of the bidomain model for some different 
values of o; and ce. In the left plot, o is fixed at the value specified in Table 1.1, and in the right 
plot ст; is fixed at the value specified in Table 1.1. 


10.3 The Monodomain Model 


As mentioned above, the bidomain model is often referred to as the gold-standard 
for computational cardiac electrophysiology. But in many cases, results of relevant 
accuracy can be achieved by using the simpler monodomain model. Here will show 
that, under one specific condition, the two models actually give the same results. We 
can show this by considering the linear PDE that needs to be solved in each time 
step; see (10.32) and (10.33) above. We assume that the conductivities are constant 
in space, and that they are related as follows, 


Ss i jay 4 J (1041) 


where J is a positive constant. By using this assumption, we note that the second 
equation of the system 


2 = У. (охур) + V - (o; Vue), (10.42) 


0 2 У. (о; Уо) + У. ((o; + Ge) Vue), (10.43) 
can be rewritten as follows, 
0 2 V-(o;Vv) + (1 + A)V - (o; Vue), (10.44) 
and therefore, 
1 
V - (mjVug) = -——~V - (a;Vv). 10.4 
(c; Vue) EF (c; Vv) (10.45) 


By inserting this observation in (10.42), we observe that we can remove ue from this 
equation and get the following scalar equation (only v is unknown here), 
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д 4 
uL Е 


аА (91%). (10.46) 


10.3.1 Operator Splitting for the Monodomain System 


In the case of the bidomain model, operator splitting involved alternating solution of 
the ODEs by the scheme (10.34)-(10.36) and the PDEs by solving the linear system 
(10.37) and (10.38). For the monodomain case, the PDE part can be simplified 
considerably by solving the system 


п+1 _ ,n41/2 А 

x6 x = An, (10.47) 
(10.48) 

which can be written to the computational form 
(I — yAtA;)u"*! = v", (10.49) 
(10.50) 

where 
Е 4 
ETIA 


10.4 Comments and Further Reading 


1. Introductions to the bidomain and monodomain models can be found in, e.g., 
[4, 16]. 

2. Numerical simulations of cardiac electrophysiology based on the bidomain and 
monodomain equations have been applied in a very large number of papers. Since 
1990, a steady stream of results have been produced by many excellent research 
groups. There are far too many studies for us to review here; instead we will refer 
you to a few papers that have become classics in the field; [3, 5, 12, 13, 18, 19]. 

3. Numerical solution of the Poisson equation is one of the best studied problems 
in scientific computing. The theory of fast linear solvers is very well developed 
for symmetric and positive definite linear systems. The linear systems arising 
from the bidomain model represent extensions of the classical systems generated 
from the Poisson equation and have therefore received substantial interest. In 
the system (10.39), the matrices A; and Ae are symmetric and thus the complete 
system can be written in a symmetric form if the lower part of the system is 


multiplied by =. . Then the system reads 
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Мод. LOAD A. п+1 n+1/2 
f уо? кез im _ P +1/ |. (4051) 
“+С, Ai “+С, (А; + Ae) и 0 


This system is symmetric and positive definite and therefore fast solvers can be 
applied; see, e.g., [2, 4, 7, 8, 10, 11, 14, 15, 16, 17]. 

4. The units used in spatially resolved electrophysiology models can be horrible 
to keep track of and it is not uncommon to miss the target by a factor of 10? or 
even 10°, and little sense arises from such computations. After some time you 
will get used to this difficulty and will at least become wise enough to hide your 
blunders until they are properly corrected. As usual, there is no way around this 
but blood, toil, sweat and tears. However, finite differences can actually be of 
some assistance. We find it easier to check the units when a system of equations 
is written in the form of a finite difference scheme, because the derivatives in 
the original equations often cause confusion, and differences are easier to deal 
with. Note also that the basic rules are quite simple; all terms in a sum must be 
expressed in terms of the same unit, and so must the left- and right-hand sides 
of an equation. 

5. There are two main reasons for approximating the bidomain model by the 
monodomain model. The first reason is the complexity of solving the equations. 
The monodomain model is of a very classical form and software are available 
from other fields and in general finite element libraries. The bidomain model, on 
the other hand, is less standard and implementing solution methods is therefore 
considered more challenging. Note, however, that excellent open-source software 
libraries are available; see, e.g., https: //opencarp.org. The second reason is 
that the computational complexity (CPU-efforts needed to solve the equations) 
is often regarded to be much higher for the bidomain model than for the the 
monodomain model. This strongly depends on the methods used to solve the 
equations; see, e.g., [17]. 

6. It has been shown that solutions of the monodomain model, in many cases, 
provide very good approximations of the solutions of the bidomain model, see, 
e.g., [12]. 

7. The bidomain model is often hailed as the gold standard of computational 
electrophysiology, and we have added to this acclaim. Forty years ago, it was 
almost unthinkable to perform simulations of a whole heart because of the 
computational complexity. In 1984, it was estimated that it would take 3000 
years to solve the bidomain model for 10 ms using a mesh with a million 
nodes [1], whereas the representation of the full human heart required about 26 
million nodes. This estimate was clearly on the pessimistic side, since in 2006, 
the full simulation was performed in only two days (see [12]), and a few years 
later such simulations could be performed in minutes (see, e.g., [9]). Today, 
the bidomain model is used routinely and simulation times are acceptable even 
without extreme computing facilities. 
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A 
Check for 


Chapter 11 се 


Re-Introducing the Cell: The 
Extracellular-Membrane-Intracellular (EMI) 
Model 


As mentioned earlier, the bidomain system is currently the standard mathematical 
model of cardiac electrophysiology. This system is now routinely solved and provides 
valuable insights into the conduction of electrical signals in cardiac tissue. However, 
the model has one glaring limitation: The cardiomyocyte is nowhere to be found in 
the model, since the extracellular space, the intracellular space and the cell membrane 
are all assumed to be everywhere in the computational domain. The cell was lost in 
homogenization! There is a tremendous advantage to this because the model becomes 
much simpler and thus solvable for the whole human heart. And it works! But the 
downside is of course that the cell is the essential building block of the tissue and 
leaving it out of the model has consequences. For instance, it becomes impossible 
to investigate the detailed dynamics of the electrochemical processes in the vicinity 
of a small collection of cells. 

Here, we will present an alternative cell-based model. The model represents the 
extracellular (E) domain, the cell membrane (M) and the intracellular (I) domain 
explicitly and it is therefore referred to as the EMI model. The main advantage of 
this model is that it becomes feasible to represent the cell and the cell membrane 
in a much more detailed manner. For instance, it is possible to study both the effect 
of varying ion channel density across the cell membrane and cell to cell variations 
of properties in the model. But it comes with a stiff price: both implementation and 
computing efforts are much more demanding than for the bidomain model. 

Here, we will present the EMI model and observe that, again, we can come 
up with a reasonable solution method by applying operator splitting and replacing 
derivatives with differences. We will present a case that is as simple as possible but 
also give references to more challenging applications. 
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Me OQ. 


Fig. 11.1 Illustration of the different parts of the domain in the EMI model. A cell, О; is surrounded 
by an extracellular space, Qe. The interface between О; and Qe defines the cell membrane, Г. Note 
that all the EMI model simulations are performed in 3D. 


11.1 The EMI Model for one Cell Surrounded by an 
Extracellular Space 


The system of partial differential equations forming the EMI model for a single cell 
surrounded by an extracellular space like illustrated in Fig. 11.1 is given by (see, 
e.g., [1, 2, 14, 41]): 


V-o;Vu; = 0, in Qj, (11.1) 

Vos Vue = 0, in О„, (11.2) 
Ue = 0, at 0, (11.3) 

De : C2 Vue = =; зо Уи, at T, (11.4) 
Hj = Ue = 0, atl, (11.5) 

Im = -nji-:o;Vui, at T, (11.6) 

2 = cs — lion), аїГ. (117) 


Here, the unknown variables to be found are the intracellular potential, u;, the 
extracellular potential, ue, and the membrane potential, v, all given in units of 
millivolts (mV). The intracellular potential is defined in the intracellular space, О;, 
the extracellular potential is defined in the extracellular space, Qe, and the membrane 
potential is defined at the membrane, Г, defined as the interface between Q; and Qe. 
The outer boundary of the extracellular space is denoted by 0€2,. Time is given in 
milliseconds (ms) and distance is specified in centimeters (cm). Furthermore, o; is 
the intracellular conductivity (in mS/cm), сг is the extracellular conductivity (in 
mS/cm), Cm is the specific membrane capacitance (in uS/cm?), and п; and ne are 
the outward pointing unit normal vectors of Q; and Qe, respectively. Like in the 
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membrane models in Chapter 8, the cable equation in Chapter 9 and in the bidomain 
and monodomain models in Chapter 10, the term Jion represents the current density 
(in uA/cm?) through ion channels on the cell membrane. These current densities 
can, for example, be modeled by the Hodgkin-Huxley model described in Chapter 8. 
In that case, 

lion = ма + Ik + h, (11.8) 


and we get the following extra equations for the state variables defined at the 
membrane, Г: 


от = ay (1 — m) – Bm, atT, (11.9) 
оп = оһ(1— h) – Bah, atT, (11.10) 
x =a,(1-r)- br, aT. (11.11) 


11.1.1 Numerical Scheme for the EMI Model 


As observed in the previous chapters, a numerical finite difference scheme for the 
EMI model can be defined by applying operator splitting and replacing derivatives 
with differences. We define a scheme where for each time step, n, there is one 
unknown, uz, for all mesh points in the extracellular space and one unknown, и”, for 
each intracellular point. In addition, for the mesh points located on the membrane, 


there are six unknowns, ий, ие, v", m", h”, andr”. 


Operator Splitting for the EMI Model 


We define an operator splitting scheme for the EMI model where for each time step 
we first solve the nonlinear ordinary differential equation part of the system 


Ov 1 

Se tT, 11.12 
Ot Cn И ( ) 
a = ay(1 — m) — бт, at Г, (11.13) 
oh 

ЕП = ол(1—Л)— fih, at Г, (11.14) 
x = а,(1- ғ) – fr, аї Г, (11.15) 


with initial conditions from the previous time step. Then, in the second step of the 
operator splitting scheme, we solve the linear system, 
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У. ои; = 0, in О;, (11.16) 

V - OeVue = 0, in Qe, (11.17) 
Ue = 0, at Qe, (11.18) 

De : Oe Vue = —nj · 07 VUj, агг, (11.19) 
Ui — Ue = U, аї Г, (11.20) 

Im = —ni: сУи, at Г, (11.21) 

2 = col аї Г, (11.22) 


with initial conditions provided from the first step of the operator splitting scheme. 


Finite Difference Approximation of the EMI Model 


The first step of the operator splitting scheme for the EMI model is simply a system 
of ordinary differential equations, and we find the numerical solutions by replacing 
the derivatives in the form д with the differences in the form _ in an explicit 
manner. 

In the second step of the operator splitting scheme, we also replace the temporal 
derivative in (11.22) with the standard difference, but we here treat the system in an 


implicit manner. That is, we replace (11.22) with 


п+1 п 1 1 

———— = Ш}, 11.23 
At Ca ( ) 
Furthermore, we use standard differences for the derivatives in (11.16), (11.17), 
(11.19), and (11.21). However, some special treatment is required for the normal 
derivatives in (11.19) and (11.21) at the corners of the cell. The details of the 
finite difference scheme is found in the code associated with these notes and is also 
described in more detail in [41]. 


11.1.2 EMI Model Simulation of a Neuronal Axon 


Using the parameter values specified in Table 11.1, we perform an EMI model 
simulation of а neuronal axon using a similar setup as for the cable equation in 
Chapter 9. However, to reduce the computational cost, we consider a shorter axon 
than in Chapter 9 (0.2 cm). A traveling wave moving from left to right is initiated 
by increasing the membrane potential of the leftmost 0.05 cm of the axon. Fig. 11.2 
shows the numerical EMI model solution of the problem at three different points in 
time along a plane in the x- and y-directions. The left panel shows the extracellular 
potential, and the right panel shows the membrane potential, v. 
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Table 11.1 Parameter values used in the EMI model simulations of an axon. The parameter values 
of the Hodgkin-Huxley model are as specified in Chapter 8. 


Parameter Value Parameter Value 
Gn 1 uF/cn? О; 2000 um x 10 um x 10 um 
Ci 4 mS/cm Qi U Qe 2060 um x 50 um x 50 um 
Te 3 mS/cm Ax 10 um 
At 0.02 ms Ay, Az 2.5 um 
T Extracellular potential Membrane potential 
0.04 
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Fig. 11.2 EMI model solution for three points in time along a plane in the x- and y-directions. 
The plane is located at the z-value corresponding to the upper boundary of the cell. The left panel 
shows the extracellular potential, and the right panel shows the membrane potential, v. 


11.2 The EMI Model for Connected Cardiomyocytes 


The EMI model for one cell considered in the previous section can be extended to 
incorporate currents between individual cells through gap junctions and thus be used 
to model collections of connected cardiomyocytes (see, e.g., [14, 16, 29, 35, 36, 37]). 
The resulting system of equations can be solved using a similar operator splitting 
technique as the one applied for a single cell above (see, e.g., [40]). Furthermore, a 
spatial operator splitting approach can be introduced in order to split the linear part 
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Fig. 11.3 EMI model solution (intracellular potential, иг) for a pulmonary vein sleeve at ten points 
in time, following the simulation set-up used in [20]. The two mutations N588K and A130V are 
both present. The solutions are found using the finite element method (see [20]). The finite element 
mesh used to represent each single cardiomyocyte in the simulation is illustrated in the lower panel 
of the figure. The full cylinder of cells seen in the upper panels contains 3930 cardiomyocytes, 
each associated with about 70 computational nodes. In addition, the mesh consists of about 44,000 
extracellular nodes. 


of the system system into one system for the extracellular space and one system for 
each individual cell (see [17, 19]). 


11.2.1 EMI Model Simulation of Cardiomyocytes in the Sleeve of a 
Pulmonary Vein 


To illustrate an example of the EMI model used for connected cardiomyocytes, 
we consider a collection of myocytes located around the sleeve of a pulmonary 
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vein. In [20], the EMI model was used to study how mutations that have been 
found to be associated with increased risk of atrial fibrillation affected properties 
of the cardiomyocytes in this region, known to be a common initiation site of atrial 
fibrillation. In Fig. 11.3, we consider a collection of cells forming a cylinder around 
the vein, using the same setup as in [20]. The mesh used to represent each single 
cell in the simulation is illustrated in the lower panel of the figure. The properties 
of the individual cells vary according to known differences for cardiomyocytes in 
this region (see [20]). In the simulation, a combination of two mutations found to 
be associated with atrial fibrillation is present. The first mutation, N588K, leads 
to an increased potassium current (/k;), and thus shortening of the action potential 
duration, whereas the second mutation, A130V, leads to reduced sodium current 
(Una) and thereby reduced conduction velocity. In Fig. 11.3, we observe that when 
the two mutations are present, the solution is a traveling wave continuously moving 
around the cylinder of cells. Such a reentrant excitation wave could be a potential 
mechanism of atrial fibrillation. 


11.3 Comments and Further Reading 


1. As mentioned above, the EMI model can be used to represent the individual 
cells of cardiac tissue and can therefore be referred to as a cell-based model, as 
opposed to the homogenized bidomain and monodomain models. Alternative 
cell-based models have also been introduced, including 1D single strand models 
(e.g., [6, 22, 30, 39, 42]), 2D sheet models (e.g., [9, 10, 11, 12, 21, 31, 32, 33, 34]), 
and 3D microdomain models (e.g., [24, 25, 26]). Differences between the EMI 
model and other cell-based models are discussed in [16]. 

2. The resolution used in monodomain and bidomain model simulations is often 
about Ax « 0.25 mm (see, e.g., [5, 43]). The volume of a cardiomyocyte has 
been reported to be around 16 pL, see [28]. Every computational block with 
volume (0.25 mm)? can thus cover almost 1000 cardiomyocytes (see [18]). This 
means that homogenization is very efficient in removing lots of details, which is 
good for computing efforts. But it also means that lots of details are lost, which 
is bad news for understanding the physics at the level of the cells. 

3. As mentioned above, the time to solve the bidomain model for a million nodes 
was estimated to 3000 years in 1984. So what is the estimate for EMI (now, in 
2023)? According to [18] the computing time for one time step is about 0.02 ms 
for each cell. For an action potential lasting for 500 ms, the total number of time 
steps is 500 x 10? when the time step is At — 0.001 ms. The computing time per 
cell for an action potential of 500 ms is therefore about 10 seconds. This means 
that we can easily deal with small collections of cells. Simulating 1000 cells 
would take less than three hours. But the human heart contains between 2 and 
3 billion cells (see [38]). The computing time for 2 billion cells for one action 
potential is about 2? x 10 seconds or 23,148 days or about 63 years. So it is not 
as bad as the bidomain model anno 1984, but it is a long wait! 
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What is the relation between the EMI model and the bidomain model? The 
bidomain model was developed long ago and can be derived in many different 
ways, but it can also be directly obtain by averaging the EMI equations over 
many cells (see [15]). 


. Derivations of the EMI model can be found in [1, 13, 14]. The EMI has been 


used to study a number of different electrophysiological phenomena, including 
applications relevant for both neuroscience (e.g., [2, 3, 4, 41]) and cardiac 
electrophysiology (e.g., [16, 20, 29, 35, 36, 37]). Furthermore, a number of 
numerical strategies for solving the equations have been proposed, including 
both finite difference and finite element schemes (see, e.g., [1, 17, 19, 23, 40]). 


. One of the simplifying assumptions underlying the EMI model presented in 


this chapter is that the effect of diffusion of ions in the intracellular and 
extracellular spaces are ignored. Such diffusion effects can be included in the 
model by including the Kirchhoff-Nernst-Planck (KNP) equations in the model, 
sometimes referred to as KNP-EMI models (see, e.g., [7, 8, 27]). 
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A 
Check for 


Chapter 12 wies, 
The Poisson-Nernst-Planck (PNP) Model 


In these notes, we have considered models of electrophysiology across several scales. 
The first was the membrane model. It assumes that the action potential is similar 
across the whole cell membrane, and the model represents the action potential as 
a function of time alone. No spatial variable is involved in the pure membrane 
models, so a length scale of these models does not make sense. As models of cardiac 
electrophysiology, we next considered the monodomain and bidomain equations. 
These models are accurate descriptions of the physics at the scale of millimeters, 
and the typical mesh resolution is about Ax ~ 0.25 mm (see, e.g., [1, 16]). The EMI 
model is cell-based and represents the physics at the micrometer scale. The typical 
mesh size for the EMI model is about 10 um see, e.g., [7, 8, 9]. 

We have moved from the homogenized millimeter scale (monodomain/bidomain) 
to the cell-based EMI model on the micrometer scale. Next, we move to the nanometer 
scale. The reason for this is that strong electrical and chemical gradients exist very 
close to the cell membrane. These gradients, referred to as the Debye layer (see, e.g., 
[6]), are only a few nanometers wide. In order to study what happens very close to the 
membrane, it is necessary to solve equations on the nanometer level, and the proper 
equations are referred to as the Poisson-Nernst-Planck (PNP) model. Once again 
we will see that even if the model is rather complex, we can use the standard tricks 
introduced above to solve the equations; operator splitting and finite differences are 
all we need (plus a little blood, toil, sweat and tears). 
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12.1 The PNP System of Partial Differential Equations 


The PNP system can be written in the form 


V -(eVd) = –р, (12.1) 
22 = У. руУс + У (ЮОк8кскУф), (12.2) 


for k = (Na*, K*, Ca?* апа СГ}. 


In this system, the variables are the electric potential, ф (in mV), and the ion 
concentrations, cy (in mM) for k = (Na*, К+, Ca?*, CI"). Furthermore, we have 
used the following definitions, 


E = 8,60, (12.3) 
P = po +F У) ась, (12.4) 
k 
©кё 
Eo 12.5 
Bk КЕТ (12.5) 


The parameters F, £o, £y, po, Dx, e, kg, T, and zy for k = (Na*, K+, Са?* and СГ} 
are defined in Table 12.1. 


Table 12.1 Parameter values used in the PNP simulations, taken from [10]. 


Parameter Description Value 

F Faraday's constant 96485.3365 C/mol 
50 Vacuum permittivity 8854 fF/m 

£l Relative permittivity, £+, in Q; and Qe 80 

Em Relative permittivity, =, in Qm, Qc 2 

Dya* Diffusion coefficient for Na* in Q; and Qe 1.33 - 10° nm?/ms 
Юк+ Diffusion coefficient for Kt in О; and Qe 1.96 - 10° nm?/ms 
Dco Diffusion coefficient for Ca?* in О; and Qe 0.71 - 10° nm?/ms 
Dcr Diffusion coefficient for CI” in О; and Qe 2.03 - 10° nm?/ms 
ZNat Valence of Na* 1 

ZK+ Valence of K* 1 

2ср+ Valence of Са2+ 2 

ZCI- Valence of СГ -1 

e Elementary charge 1.60217662 - 107? C 
kp Boltzmann constant 1.380649 - 107” mJ/K 
T Temperature 310K 

Ax, Ay Spatial discretization parameter 0.5 nm 

At Numerical time step 0.02 ns 


Again, we will solve this system by applying operator splitting and replacing 
derivatives by differences. We assume that the solution is known at time t, = nAt. The 
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Table 12.2 Initial conditions for the ion concentrations in the intracellular space, Q;, and the 
extracellular space, Qe. In the membrane, Q,,,, all ion concentrations are set to zero. Furthermore, 
in the K* channel embedded in the membrane, the concentration of K* ions varies linearly from 
the intracellular to the extracellular concentration, whereas the remaining ion concentrations are 
set to zero. To make the entire domain electroneurtral at t = 0, po is set up so that p = 0 in the 
channel for the initial conditions (see (12.4)). 


Ion Intracellular Extracellular 
Na* 12 mM 100 mM 

K* 125 mM 5mM 

Ca? 0.0001 mM 1.4 mM 

CI 137.0002 mM 107.8 mM 


first step in the algorithm is to compute the electrical potential at time „+ = (n-- 1)At. 
This step can be written as 


Va: (Vy 9*1) = — p", (12.6) 


where Ул denotes a finite difference approximation of the gradient. Note that o" is 
taken from the previous time step so only the electrical potential is unknown in this 
step. In the first time step, we use the initial conditions to compute p given by (12.4). 

When $"*! has been computed, we can compute V; ф"*! and use it to solve the 
concentration equations by the following scheme, 


c _ „п 

у = Vac Ривс + 9: (Diper! vant) (12.7) 
for k = (Na*, K+, Ca?*,CI- }. Writing the complete finite difference schemes for 
these equations is a bit messy, but the interested reader can consult the online Matlab 


code, or the supplementary information of [10]. 


12.1.1 Numerical Simulation of the Resting State 


We will use the scheme given by (12.6) and (12.7) to compute the resting state close 
to the cell membrane. In the models introduced in the previous chapters, we have 
taken the concentrations to be constants in space; i.e., we have assumed that the 
concentrations can vary in time across the cell membrane, but not in space in the 
intra- or extracellular spaces. This is often an accurate approximation, but we will see 
that significant gradients exist very close to the cell membrane. Furthermore, in the 
EMI model, the cell membrane is assumed to be infinitely thin, and in the bidomain 
and monodomain models, the cell membrane is assumed to be everywhere! In the 
PNP model simulation, the cell membrane is explicitly represented in the model 
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Fig. 12.1 Illustration of 
the PNP model domain, 
consisting of an intracellular 
domain, Q;, an extracellular 
domain, O,, a membrane 
domain, Q,,, and a Kt 
channel domain, Qo. In 
the simulation reported in 
Fig. 12.2, we use the domain 
size L; = Le = Ly = 50 nm, 
Lm =5nm, We = 5 nm. 


Fig. 12.2 The membrane 
potential, v = ф; — de, as а 
function of time in the PNP 
simulation. 


(5 nm wide), but still the ion channels integrated in the membrane are modeled in 
an oversimplified manner. 

We solve the PNP system in two spatial dimensions, see Fig. 12.1. Note that the 
computational domain consists of an intracellular domain, О;, a cell membrane, Qm, 
an extracellular space, Qe, and a К+ channel, Q,. The initial conditions are given in 
Table 12.2, and the solutions are presented in Figs. 12.2-12.4. 

In Fig. 12.2, we have plotted the membrane potential (v = ф; — фе) as a function 
of time during the PNP simulation, and we see that the value starts at 0 and gradually 
approaches a typical resting potential value of about —80 mV. In Fig. 12.3 and 
Fig. 12.4, we show how the PNP model solution varies in space at the end of the 
simulation (at t = 50 ns) in the part of the domain that is located in the 5 nm closest 
to the membrane. The upper panel focuses on the intracellular side and the lower 
panel focuses on the extracellular side of the membrane. In Fig. 12.3, we plot the 
full 2D solution, and in Fig. 12.4, we plot the solution along lines in the x-direction 
at y = 0 nm and y = 25 nm. We observe that a boundary layer is formed close to the 
membrane with slightly different values of the ion concentrations and potential than 
in the bulk intracellular and extracellular spaces. 
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Fig. 12.3 Solutions at the end of the simulation (t = 50 ns) of the PNP model simulation in the 
5 nm closest to the membrane on the intracellular side (upper panel) and the extracellular side 
(lower panel). 


Ф + + - 
ф Ма K „лоз Ca? cl 
-82.8 ————_, 125 138.5 
== yeon 12 10 
-82.9 у= 290 9.95 138 
=> > > 124.5 > > 
ct E 11.95 E E 99 E 
-83 137.5 
m 9.85 
-83.1 11.9 98 137 
46 48 46 48 46 48 46 48 46 48 
0 101.5 54 144 108 
101 53 1.43 107.5 
Ф > 02 = = = 142 = 
Сі Е Е Е Ё 107 
100.5 52 44i 
106.5 
-04 14 
100 54 466 
56 58 60 56 58 60 56 58 60 56 58 60 56 58 60 
х (пт) х (пт) х (пт) х (пт) x (nm) 


Fig. 12.4 Solutions at the end of the simulation (t = 50 ns) of the PNP model simulation in the 
5 nm closest to the membrane on the intracellular side (upper panel) and the extracellular side 
(lower panel). We show the solutions along lines in the x-direction for y = 0 nm and y = 25 nm. 
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1. The PNP equations modeling the electrical potential and ionic concentrations 
close to biological membranes have been studied by several authors; see, e.g., 
[4, 5, 10, 11, 14]. But the PNP equations are also used to model lithium ion 
batteries, see, e.g., [12, 17]. 

2. The model, methods and setup in this chapter was motived by the paper [10]. 

3. A simplified version of the PNP equations are referred to as the KNP 
(Kirchhoff-Nernst-Planck) equations. In these equations, electroneutrality is 
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assumed everywhere, meaning that p = О in the entire domain. That is, 
the charges sum to zero everywhere. Numerical approximations of the KNP 
equations can be found in, e.g., [2, 3, 13, 15]. 


. Above, we mentioned that one single computational block (Ax? = (0.25 mm)?) 


for a standard mesh used to solve the bidomain model is large enough to 
cover almost 1000 cardiomyocytes. In the PNP model we use the resolution 
Ax = 0.5 nm and thus the volume of one block is 0.125 nm?. The volume 
of a sodium atom is ~ 0.0244 nm? so one computational block covers about 
five sodium atoms. The next scale, following bidomain/EMI/PNP, is therefore 
simulation based on representation of individual atoms. If a cell with a volume 
of 16 pL is represented by a uniform mesh at atomic (sodium atom) resolution, it 
will require about 6.5 x 10! blocks, which is a lot! A reasonably well-equipped 
PC today has 16 GB memory and can therefore work with a vector (in Matlab) 
of ~ 2 x 10? real numbers. Thus, about 325,000 of these PCs are needed to store 
one real number per atom in a cardiomyocyte of 16 pL. So, it will probably take 
some time before atomic scale simulations can be used to simulate whole cells. 
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